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SUMMARY 


In  this  report  the  effect  of  multiple  scattering  on  the  coherent  wave 
propatation  in  discrete  random  media  has  been  investigated.  The  media  is 
modelled  as  a  random  distribution  of  spherical  and  non-spherical  scatterers. 
The  first  and  second  order  probability  distribution  functions  are  specified 
and  a  self-consistent  T-matrix  approach  together  with  Lax's  guasi- 
crystalline  approximation  is  used  to  derive  dispersion  equations  whose 
singular  solutions  yield  the  complex  propagation  constant  of  the  "effective" 
medium. 

Our  formalism  is  well  suited  for  numerical  computations  for  wavelengths 
comparable  to  scatter  size  and  high  volume  concentration  of  scatterers;  to 
our  knowledge,  it  is  the  only  method  that  provides  reliable  numerical 
results  for  the  attenuation  and  phase  velocity  of  the  coherent  wave.  The 
quasi-crystalline  approximation  (QCA)  is  found  to  be  applicable  for  all 
concentration  values  from  low  to  high.  The  QCA  that  is  used  to  truncate  the 
hierarchy  of  equations  during  the  configurational  averaging  procedure 
requires  knowledge  of  the  two  particle  joint  probability  distribution 
function. 

Our  model  of  the  random  system  is  that  the  sphere  circumscribing  the 
scatterers  cannot  interpenetrate.  In  the  statistical  mechanics  literature 
this  is  synonymous  with  ensemble  of  'hard  spheres . '  For  such  a  model, 
several  forms  of  the  pair-correlation  function  can  be  obtained  from  the 
statistical  mechanics  of  simple  liquids  by  using  several  theories  and 
calculations  such  as  the  Hypernetted-chain  Equation  (HNC) ,  the  Percus- 
Yevick  Approximation  (P-YA),  the  Self-Consistent  Approximation  (SCA), 

Monte  Carlo  calculations,  etc.  For  the  hard  sphere  model,  temperature 
and  other  thermodynamic  quantities  do  not  appear  in  the  final  form  of  the 
correlation  functions  and  hence  are  equally  valid  for  our  system.  Five 
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forms  of  the  pair  correlation  function  were  considered: 

(a)  Well  stirred  approximation  (WSA) .  This  is  the  simplest  and  in¬ 
corporates  only  the  whole  correction.  In  the  Rayleigh  limit,  for 

c  >  *125,  this  model  fails.  But  with  increasing  frequency,  it  is  good 
for  higher  values  of  c.  At  ka  =  ua/c  -  5.0,  there  is  no  difference 
between  this  and  other  models  suggesting  that  correlations  are  unimportant 
at  high  values  of  the  wavenumber. 

(b)  The  Percus-Yevick  approximation  can  be  used  to  obtain  a  semi- 
analytical  form  of  the  pair  correlation  function.  It  is  good  for 

c  <  0.35. 

(c)  Virial  Expansion  in  powers  of  the  number  density  can  also  be 
conveniently  used  since  analytical  expansions  are  available.  Obviously 
it  works  only  for  c  <  *1. 

(d)  The  Matern  model  is  an  analytical  model  for  the  pair  correlation 
function  that  is  also  convenient  to  use  for  c  <  *125.  The  Matern  model  is 
found  to  be  superior  to  the  WSA  for  c  <  0.125. 

(e)  The  self-consistent  approximation  to  the  pair  correlation 
function  which  is  based  on  the  Percus-Yevick  model  and  the  Hyper-netted 
Chain  approximation  is  the  one  that  has  proved  most  successful  in  our 
computations.  It  provides  reliable  results  for  a  wide  range  of  scatterer 
concentrations . 

All  of  the  above  forms  were  tested  for  several  values  of  c  and  ka. 
Recently  our  calculations  were  compared  with  the  experimental  findings  of 
Professor  A.  Ishimaru  for  a  distribution  of  latex  spheres  in  water.  The 
agreement  is  excellent  for  all  cases,  see  our  paper  Nos:  1,  3  and  4.  In 
paper  Nos:  4  and  5,  we  have  introduced  the  concept  of  an  average  T- 
matrix  to  include  a  size  distribution  for  scatterers  using  the  Gaussian 
distribution  function.  The  question  may  be  raised  as  to  the  radius  of  the 
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excluded  volume  surrounding  each  scatterer.  If  we  use  4a,  the  diameter  of 
the  largest  size  sphere  of  the  distribution,  although  we  would  be  correct  in 
the  integration  on  the  allowed  position  for  each  scatterer,  we  would  be 
excluding  an  unreasonable  volume,  thus  limiting  the  volume  fraction  of 
scatterers.  It  seems  more  reasonable  to  take  the  diameter  of  the  mean 
sphere,  2a  as  the  radius  of  the  excluded  spherical  volume.  On  the  average, 
this  would  be  applicable  to  most  scatterers  and  allow  us  to  consider  higher 
volume  fractions.  If  the  volume  fraction  is  low  enough,  whether  the  radius 
of  the  excluded  sphere  is  2a  or  4a  will  not  pose  a  problem,  see  our  paper 
No:  6  for  more  details.  In  Paper  No.  5,  we  have  compared  Keller's  and 
Twersky's  multiple  scattering  approaches  with  our  formalism.  Our  formalism 
is  in  exact  agreement  with  Twersky's  approach.  However,  it  is  better 
suited  for  numerical  computations.  Recently,  we  have  also  performed  cal¬ 
culations  on  coated  dielectric  spheres  in  free  space  or  in  another 
dielectric.  Numerical  results  are  obtained  for  the  complex  average 
dielectric  constant,  coherent  attenuation  and  phase  velocity  as  a  function 
of  concentration  (0£c<^0.42).  A  plot  of  coherent  attenuation  vs  c  for 
coated  spheres  is  attached  with  this  report,  see  Fig.  1.  A  T-matrix  program 
has  also  been  developed  to  obtain  the  T-matrix  of  coated  irregular  arbitrary 
shaped  bodies  which  can  be  implemented  in  our  multiple  scattering  approach 
to  study  wave  propagation  through  snow,  aerosol  and  various  kinds  of  debris. 
Numerical  results  of  irregular  arbitrary  shaped  bodies  will  be  forthcoming. 

Recently,  we  had  extended  our  formalism  to  study  high  frequency 
propagation  of  waves  in  random  media.  We  have  compared  our  results  with 
those  of  experimental  findings  obtained  by  Ishimaru  for  ka  as  high  as 
83.596  and  the  agreement  is  found  to  be  excellent,  see  Fig.  2.  We  are 
confident  that  we  have  a  well  founded  formalism  and  a  sound  numerical 
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technique  along  with  our  T-matrix  and  finite  element  (uni-moment)  approaches 
for  the  study  of  wave  propagation  in  random  media. 

Various  publications  resulting  from  our  investigations  are  enclosed. 
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Fig.  1.  Attenuation  Im(K/k^)  as  a  function  of  concentration  c  of 
/“  coated  spheres  (  a./a-  =  1.17  )  for  three  values  of  the  free  space 

wavenumber  kQa|  2 
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Fig.  2.  Plot  of  Y  =  2  Im  K/n^  versus  concentration  where  nQ 
is  the  number  density  and  is  the  extinction  cross  section  of 
a  single  sphere  with  ka  =  83.596. 
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Coherent  electromagnetic  wave  propagation  in  an  infinite  medium  composed  of  a  random  distri¬ 
bution  of  identical,  finite  scattercrs  is  studied.  A  self-consistent  multiple  scattering  theory  using  the  T 
matrix  of  a  single  scatterer  and  a  suitable  averaging  technique  is  employed.  The  statistical  nature  of  the 
position  of  scalterers  is  accounted  for  by  ensemble  averaging.  This  results  in  a  hierarchy  of  equations 
relating  the  different  orders  of  correlations  between  the  scalterers.  Lax's  quasi-crystalline  approxi¬ 
mation  is  used  to  truncate  the  hierachy  enabling  passage  to  a  homogeneous  continuum  whose  bulk 
propagation  characteristics  such  as  phase  velocity  and  coherent  wave  attenuation  can  then  be  studied. 
Three  models  for  the  pair  correlation  function  are  considered.  The  Matem  model  and  the  well-stirred 
approximation  are  good  only  for  sparse  concentrations,  while  the  Percus-Yevick  approximation  is 
good  for  a  wider  range  of  concentration.  The  results  obtained  using  these  models  are  compared  with 
the  available  experimental  results  for  dielectric  scatterers  embedded  in  a  host  dielectric  medium.  Practi¬ 
cal  applications  of  this  study  include  artificial  dielectric  (composites)  and  electromagnetic  wave  propa¬ 
gation  through  hydrometeors.  dust,  vegetation,  etc. 


1.  INTRODUCTION 

We  consider  the  propagation  of  plane  coherent 
Q|  electromagnetic  waves  in  an  infinite  medium  contain- 
•»■-  ing  identical,  lossless,  randomly  distributed  particles. 

Our  aim  is  to  characterize  the  random  medium  by  an 
..  effective  complex  wave  number  K  which  would  be  a 
v"\  function  of  the  particle  concentration,  the  electrical 
size,  and  the  statistical  description  of  the  random 
positions  of  the  scatterers.  The  imaginary  part  of  K 
describes  the  coherent  attenuation  which  is  due  to 
multiple  scattering  only  when  the  particles  them¬ 
selves  are  assumed  to  be  lossless.  The  understanding 
of  the  behavior  of  Im  (K)  as  a  function  of  particle 
■jj  concentration  (c)  and/or  frequency  (ka)  is  very  im¬ 
portant  in  many  practical  applications,  including 


].’•  'Now  at  Colorado  State  University,  Fort  Collins,  Colorado 
80523. 

Copyright  1982  by  the  American  Geophysical  Union. 

kl 

Paper  number  2S04I I. 

.  0048-6604  82  0910-041 1  $08  00 


wave  propagation  in  the  atmosphere  and  oceans  and 
whenever  random  distributions  of  scatterers  influence 
electromagnetic  wave  behavior. 

The  theoretical  formulation  presented  here  closely 
follows  the  procedure  described  by  Varadan  et  al. 
[1979,  1982],  Varadan  and  Varadan  [1980],  and 
Varadan  [1980].  This  approach  is  based  on  a  self- 
consistent  multiple  scattering  theory  and  relies  on  the 
T  matrix  [ Waterman ,  1971]  which  relates  the  field 
scattered  by  a  particle  to  an  arbitrary  exciting  field. 
The  statistical  description  of  the  random  positions  of 
the  scatterers  is  used  to  define  a  configurational 
average  which  results  in  a  hierarchy  of  equations 
relating  the  different  orders  of  correlations  between 
the  scatterers.  Lax's  [1952]  quasi-crystalline  approxi¬ 
mation  (QCA)  is  used  to  truncate  the  hierarchy 
which  results  in  the  usual  ‘hole  correction’  integrals. 
Following  Twersky  [1977,  1978a.  fr],  a  radially  sym¬ 
metric  pair  correlation  function  is  introduced,  and 
approximate  models  are  chosen  from  Talbot  and 
Willis  [1980].  The  ‘well-stirred’  approximation 
(WSA)  which  was  used  previously  by  Varadan  et  al. 
[1979]  and  Bringi  et  al.  [1981]  assumes  no  corre- 
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Fig.  I.  Geometry  of  randomly  distributed  and  aligned  scatterers. 


lation  between  the  particles  except  that  they  should 
not  interpenetrate.  In  particular,  the  WSA  gives 
unphysical  results  for  0.125  at  the  Rayleigh  or 
low-frequency  limit. 

In  this  paper  we  consider  two  other  pair  corre¬ 
lation  functions,  viz.,  (1)  the  Matern  [1960]  model 
and  (2)  the  Percus-Yevjck  [ Percus  and  Yeviek,  1958] 
model  for  a  classical  system  of  hard  spheres.  Compu¬ 
tations  of  Im  ( K )  are  presented  for  dielectric  scat¬ 
terers  imbedded  in  a  host  dielectric  medium,  using 
the  above  three  models  as  a  function  of  frequency 
and  concentration.  We  compare  our  solutions  to 
some  recent  optical  propagation  experiments  con¬ 
ducted  by  Ishimaru  (A.  Ishimaru,  personal  communi¬ 
cation,  1981).  Sample  computations  are  also  presen¬ 
ted  comparing  the  WSA  and  the  single  scattering 
approximation  for  a  rain  medium. 

2.  FORMULATION  OF  THE  PROBLEM 

Consider  (V  identical,  finite  dielectric  scatterers 
that  are  randomly  distributed  either  in  free  space  or 
in  a  host  dielectric  medium.  The  scatterers  are  homo¬ 
geneous  with  a  relative  dielectric  constant  of  e,,  their 
centers  being  denoted  by  0,,  02,  03,  ■  •  • ,  0„  (see 
Figure  I).  They  are  assumed  to  be  bodies  of  revol¬ 
ution  with  symmetry  axis  parallel  to  the  z  direction. 
A  monochromatic,  plane,  coherent  electromagnetic 
wave  is  assumed  to  propagate  along  the  symmetry 
axis  of  the  scatterers  so  as  to  satisfy  the  condition 
that  the  effective  medium  be  isotropic  and  polariza¬ 
tion  insensitive.  The  time  dependence  of  the  incident 
field  and  hence  the  fields  scattered  by  the  individual 


scatterers  is  of  the  form  exp  ( — jiot )  and  is  suppressed 
in  the  equations  that  follow.  Even  though  the  theory 
presented  here  is  valid  for  spheroidal  scatterers 
[  Varadan  et  ai,  1982],  we  present  computations  only 
for  spherical  scatterers  in  order  to  compare  our  re¬ 
sults  with  available  experiments. 

Let  E°(r)  be  the  electric  field  arising  from  the  inci¬ 
dent  plane  wave  and  E?(r)  the  field  scattered  by  the 
ith  scatterer.  Both  these  fields  satisfy  the  vector 
Helmholtz  equation.  The  total  field  at  any  point  out¬ 
side  the  scatterers  is  given  by  the  sum  of  the  incident 
field  and  the  fields  scattered  by  all  the  scatterers, 
which  can  be  written  as 

S 

E(r)  =  E°(r)  +  £  E?(p,)  Pi  =  r  -  r(  (1) 

i  =  1 

where  Ef(pf)  is  the  field  scattered  by  the  ith  scatterer 
at  the  observation  point  r.  However,  the  field  that 
excites  the  ith  scatterer  is  the  incident  field  E°  plus 
the  fields  scattered  from  all  other  scatterers  except 
the  ith.  The  term  exciting  field  E'  is  used  to  dis¬ 
tinguish  between  the  field  actually  incident  on  a  scat¬ 
terer  and  the  external  incident  E°  produced  by  a 
source  at  infinity.  Thus  at  a  point  r  in  the  vicinity  of 
the  ith  scatterer,  we  write 

.V 

E'(r>  =  E°(r)  +  £  Ej  <p;)  a  <;  |  Pj  \  <  2a  (2) 

i*i 

where  a  is  a  typical  dimension  of  the  scatterer. 

The  exciting  and  scattered  fields  for  each  scatterer 
can  be  expanded  in  terms  of  vector  spherical  func- 
tionswith  respect  to  an  origin  at  the  center  of  that 
scatterer: 

Ef(r)  =  X  £  £  X  b\lm„  Re  \|»/M,<p() 

r  -  1  1*1  m  *  0  a  =  i 

=  I  b'„  Re  v|»;„  (3) 

Ml 

and 

E,s(r)  =  £  (4) 

Ml 

where  n  represents  a  combined  index  notation  for  the 
set  [i.m.ff]. 

The  vector  spherical  functions  are  defined  as 

'I'n.jr)  =  V  x  [r/»,(fcr)']V',1M(0.  <p)  (5) 

=  (1  MV  x 
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In  (3-6),  k  is  the  wave  number,  h,  is  the  spherical 
Hankel  function  of  the  first  kind,  and  the  Ylma(B ,  <p) 
are  the  normalized  spherical  harmonics  defined  with 
real  angular  functions.  In  (3)  the  exciting  field  is  ex¬ 
panded  in  terms  of  the  regular  (Re)  basis  set  (Re 
obtained  by  replacing  h,  in  (S)  and  (6)  by  jH,  the 
spherical  Bessel  functions  of  the  first  kind.  Thus  the 
choice  of  the  basis  set  in  (4)  satisfies  the  radiation 
condition  at  infinity  for  the  scattered  field,  while  the 
choice  in  (3)  satisfies  the  regular  behavior  of  the  ex¬ 
citing  field  in  the  region  a  <  |p{|  <  2a.  The  super¬ 
script  i  on  the  basis  functions  refer  to  expansions 
with  respect  to  0,,  and  i>}„  and  B'„  are  the  unknown 
exciting  and  scattered  field  coefficients,  respectively. 
We  also  expand  the  incident  field  in  terms  of  vector 
spherical  functions: 

E°(r)  =  e‘u  "  £  a„  Re  «kUPi)  (7) 

m 

where  the  a,„  are  the  known  incident  field  coef¬ 
ficients. 

The  unknown  coefficients  bj.  can  be  related  to  B}„ 
by  means  of  any  convenient  scattering  operator;  in 
this  case  we  employ  the  T  matrix  as  defined  by 
Waterman  [1971]: 


IT}., 


u 

tw  ur>n' 


the  following  set  of  algebraic  equations  for  the  excit¬ 
ing  field  coefficients  b‘,„ : 

V 

b'„  =  e‘ki  r,a„  +  £  I  I  <W  JP,,> 

j  4  i  f'd»  z»n» 

•  (in 

From  (11)  it  can  be  seen  that  the  exciting  field 
coefficients  of  the  ith  scatterer  explicitly  depend  on 
the  position  and  orientation  of  the  other  scatterers. 
In  this  paper  we  consider  a  random  distribution  of 
spherical  scatterers  and  the  case  when  both  x 
and  the  total  volume  accessible  to  the  scatterers 
V  — *  x  such  that  N/V  =  n0  is  a  finite  number  den¬ 
sity.  For  such  distributions  a  configurational  average 
of  (11)  can  be  made  over  the  positions  of  all  the 
scatterers  [see  Varadan  et  al.,  1982],  The  quasi- 
crystalline  approximation  [Lax,  1952]  can  then  be 
invoked  to  arrive  at  an  equation  for  the  configur¬ 
ational  average  <6}„>i  of  the  exciting  field  coefficients 
with  one  scatterer  fixed: 

<*>!„>.  =  r‘a„  +(N-  1)  I  £ 


x r  MjKbi >,  dr, 


Substituting  (3),  (4),  and  (7)  in  (2),  we  obtain 
£  6}.  Re  <K.  =  £  a,.  Re  ♦}.+  ££  B{.  (9) 

tn  rn  ]  4  f  m 

Since  the  field  quantities  are  expanded  with  respect 
to  centers  of  each  scatterer,  we  obtain  (9)  with  basis 
functions  expanded  with  respect  to  ith  and  yth  cen¬ 
ters.  In  order  to  express  them  with  respect  to  a 
common  origin  0.,  we  employ  the  translation  addi¬ 
tion  theorems  for  the  vector  spherical  functions  (see, 
for  example,  Bostrom,  [1980])  which  can  be  written 
in  a  compact  form  as  follows: 

'K.(pj)  =  I  v„.<P,j)  Re  ♦,.,.(Pi)  I  Pul  >  1  Pi  1 

tw  (10) 

*  Jp j)  =  I  ,...(pW)  ■  tvJPi)  I  Py  I  <  I  Pi  I 

t'M* 

where  p,y  =  r,  -  r,  is  the  vector  connecting  0 j  to  0, , 
ar.  is  the  translation  matrix  for  the  vector  func¬ 
tions,  and  Rt,  is  a  matrix  with  spherical  Hankel 
functions  in  replaced  by  spherical  Bessel  func¬ 
tions. 

Employing  (8)  and  (10)  in  (9)  and  using  the  orthog¬ 
onality  of  the  vector  spherical  functions,  we  obtain 


where  p(r,|r{)  is  the  two-particfe  conditional  prob¬ 
ability  density. 

We  now  assume  that  the  average  or  coherent  field 
propagates  in  a  medium  with  an  effective  complex 
wave  number  K  =  (Kt  +  iK2)z  in  the  direction  of  the 
original  incident  field : 

<Ef(r)>j  =  A  exp  UK:  •  r)  (13) 

where  A  is  the  amplitude  of  the  coherent  wave.  Thus 
the  average  exciting  field  coefficients  may  be  ex¬ 
pressed  as 

<6}i™,>,  =  «*P  UK:  •  r,)yrtM,<5„,[<5tl(5»2  +  <5t2<$„]  (14) 

where  the  Kronecker  deltas  in  (14)  indicate  that  only 
the  azimuthal  index  m  —  1  contributes  while  those  in 
square  brackets  indicate  that  the  effective  medium  is 
isotropic. 

Equation  (14)  is  now  substituted  in  (12),  and  the 
conditional  probability  p(ry|r,)  is  expressed  in  the 
form 


T<r>|r->  =  i7rirr7)s/<x>  r>1 

(15) 

p(r,|r,)  =  0  x  <  1 

where  x  -  1 2a.  The  function  </(  x)  is  termed  the 
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Fig.  2.  The  Percus-Yevick  pair  correlation  function  for  hard 
spheres. 


pair  correlation  function,  which  is  assumed  to  be 
spherically  symmetric.  Equation  (12)  can  now  be  sim¬ 
plified  by  making  use  of  certain  symmetry  properties 
of  the  T  matrix  [ Waterman ,  1971]  and  by  invoking 
the  extinction  theorem  [T wersky,  1977]  to  yield  the 
following  set  of  equations  for  the  unknown  coef¬ 
ficients  Ytlma : 

-  'v — *  y  y  t 

y  t 


(  X  KK,  k,  c, 

U-M-M  J 


(16) 


where 


KK ,  *,  c,  ;.j 
6c 


[2kaj  i(2Ka)h\(2ka)  -  2Kahi(2ku)j\(2Ka)] 

(17) 


(fca)2  —  (/Ca)2 
+  24c  f  cJ[^(.x)-  l]hi(2kax)/i(2Ka.x)  dx 


L  21I/+1) 


rnr  +  in1  Tf  /  rurr  / 

‘  |_/</  + nj  Li  -i  oJlLo  o  oj 

•  [/'(/•  +  !>  +  /(/  +  i +  i  )]£„,  s„, 

+  l0  ■  lJ{[;.2-(/'-/)2]|T  + 1  + 

•  (1  —  —  (5tfl)  (18) 

In  (17),  c  =  |  na^rtQ  is  the  fractional  concentration  by 
volume,  and  in  (18), 

pi  Jz  ji  ] 

[m,  «2 

is  the  Wigner  3 -j  symbol.  If  the  integral  /  in  (17)  can 
be  evaluated  for  suitable  models  of  the  pair  corre¬ 
lation  function,  then  (16)  constitutes  a  set  of  homoge¬ 
neous,  algebraic  equations  with  Y(lma  as  the  un¬ 
knowns.  For  a  nontrivial  solution,  the  determinant  of 
the  coefficient  matrix  must  vanish.  This  yields  the 
dispersion  equation  for  the  effective,  complex  wave 
number  K.  In  the  Rayleigh  or  low-frequency  limit, 
analytical  solutions  for  K  can  be  found  while  at 
higher  frequencies  the  solutions  can  be  generated 
computationally.  The  real  part  Kl  is  the  phase  con¬ 
stant  while  the  imaginary  part  is  the  attenuation  con¬ 
stant  of  the  effective  medium. 

If  the  concentration  is  sparse,  c  -*  0,  then  g(x)  as¬ 
sumes  the  form  of  a  step  function  at  x  =  1.  This  be¬ 
havior  of  crt.x)  is  termed  the  well-stirred  approxi- 


Fig.  3.  The  coherent  attenuation  constant  ;  versus  ku  for  t,  = 
(;,(/.)  using  the  WSA. 
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mation  (WSA)  by  Talbot  and  Willis  [1980]  and  has 
been  used  previously  by  Fikioris  and  Waterman 
[1964],  Talbot  and  Willis  [1980]  also  discuss  a 
number  of  other  pair  correlation  functions  suitable 
for  use  at  higher  concentrations,  viz.,  a  model  due  to 
Matern  [1960]  and  the  model  obtained  by  solution 
of  the  Percus-Yevick  [Percus  and  Yevick,  1958]  inte¬ 
gral  equation  for  a  classical  system  of  hard  spheres 
[  Wertheim,  1963].  Other  forms  for  g(x)  can  be  used 
based  on  the  statistical  mechanical  theory  applied 
toward  the  study  of  dense  gases  and  liquids 
[Twersky,  19786],  In  this  study,  however,  we  show 
sample  computations  using  the  g{x)  corresponding  to 
the  WSA,  the  Matern  model  which  is  analytic,  and 
the  Percus-Yevick  approximation  (P-YA)  based  on 
tabulated  values  of  g(x)  given  by  Throop  and  Bear- 
man  [1965].  The  behavior  of  g(x)  versus  x  is  shown 
in  Figure  2  for  various  values  of  concentration. 

3.  COMPUTATIONS 


c 


Fig.  5.  The  normalized  attenuation  constant  y  versus  con¬ 
centration  c  for  different  values  of  ka  using  the  WSA. 


In  many  practical  applications  it  is  of  interest  to 
determine  the  concentration  levels  above  which 
multiple  scattering  effects  must  be  accounted  for  in  a 
rigorous  manner.  At  very  low  concentrations,  c  -» 0, 
the  particles  can  be  considered  as  essentially  decor- 
related  so  that  the  so-called  ‘single  scattering’  ap¬ 
proximation  (SSA)  can  be  used  leading  to  the  follow¬ 
ing  expression  for  Im  K: 

Im  (K/k)  *  $cQtJ(ka)  (19) 


where  Qeu  is  the  normalized  (with  respect  to  rta1) 
extinction  cross  section  of  a  sphere  of  radius  a. 
Brown  [1980]  has  shown  recently  that  under  the  con¬ 
ditions  of  c  — ►  0  and  no  correlation  between  particles, 
the  Foldy-Twersky  integral  equation  [see  lshimaru, 
1978]  for  the  coherent  field  reduces  to  (19)  above  and 
thus  does  not  account  for  the  effects  of  multiple  scat¬ 
tering.  To  demonstrate  that  (19)  and  the  theoretical 
procedure  given  in  this  paper  lead  to  identical  results, 
we  have  compared  computations  for  a  monodisperse 
rain  medium.  The  rain  medium  is  assumed  to  consist 
of  spherical  water  drops  at  concentrations  in  the 
range  10~4<c<10"2  which  encompasses  even  the 
heaviest  rainfall  conditions.  In  Figure  3  the  attenu¬ 
ation  constant  •/  =  4nK2;Kl  is  shown  as  a  function 
of  freqency  or  ka  using  the  WSA.  which  is  to  be 
compared  to  Figure  4  which  uses  (19)  or  the  SSA. 
The  refractive  index  of  water  is  taken  from  Ray 
[1972]  assuming  a  temperature  of  5  C.  Note  that 
both  solutions  yield  nearly  identical  results.  In  Figure 
5  we  show  •/  versus  concentration  (or  equivalently, 
rainfall  rate)  at  a  number  of  fixed  frequencies  or  ka 
values  using  the  WSA.  The  straight  line  relationship 
reflects  an  attenuation  versus  concentration  relation¬ 
ship  of  a  power  law  form.  Olsen  et  al.  [1978]  have 
used  (19)  to  derive  the  power  law  relation  between 
attenuation  and  rainfall  rate  (see  their  Figure  3). 

At  higher  concentrations,  pair  correlation  effects 
become  important.  In  the  Rayleigh  or  low-frequency 
limit.  Twersky  [19786]  has  given  an  expression  for 
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Fig.  6.  The  coherent  attenuation  Im  l  AC  k)  versus  concentration  c 
at  ka  =  0.05  for  latex  spheres  in  water. 

Im  (K,  k)  by  considering  the  leading  effects  of  the  pair 
correlation: 

Im  (AC  k)  =  1)  (£,  +  2|]3IT  (20) 

where  W  is  the  packing  factor  given  by 

IT  =  ,  =  1  +  24c  f  v:[s/(.v|-  1]  J.x  (21) 

( 1  +  2c)J  Jo 

The  above  integral  can  be  intergrated  analytically  to 
give  IT  when  t/(.x)  is  given  by  the  P-YA.  In  Figure  6 
we  show  Im  (X  k)  versus  concentration  for  a  fixed 
frequency  or  ka  =  0.05.  The  random  medium  is  as¬ 
sumed  to  consist  of  latex  spheres  (£,  =  2.26)  imbed¬ 
ded  in  water  and  corresponds  to  the  experimental 
setup  used  by  Ishimaru  (A.  Ishimaru,  personal  com¬ 
munication.  1981)  for  his  optical  propagation  experi¬ 
ments.  The  agreement  between  (20)  and  the  compu¬ 
tations  using  the  P-YA  is  excellent  as  expected  while 
both  the  WSA  and  the  Matern  model  fail  (Im  K  <  0) 
for  c  >  0.125.  However,  the  Matern  model  is  superior 
to  the  WSA  as  expected.  In  Figure  7  we  show  similar 
results  at  ka  =  0.56  for  which  Ishimaru  (personal 
communication,  1981)  has  experimental  results  (latex 
sphere  diameter  =  0.107/t,  /.  =  0.6p)  at  concentration 
values  of  0.01  and  0.1.  Note  that  the  SSA  overesti¬ 
mates  the  attenuation  constant  even  for  c  as  low  as 
0.01;  as  c  increases,  the  deviation  from  the  P-YA  is 
significant.  For  low  values  of  c.  the  P-YA.  WSA.  and 
Matern  model  are  all  in  excellent  agreement.  Even  at 
c  =  0.10.  all  three  pair  correlation  models  are  in  rea¬ 


sonably  good  agreement  with  the  experimental  value. 
For  c  >0.125  the  Matern  model,  as  well  as  the  WSA, 
deviates  considerably  from  the  P-YA.  Even  though 
experimental  values  for  c  >  0. 1  are  as  yet  unavail¬ 
able,  we  feel  the  P-YA  will  continue  to  predict  the 
correct  behavior  for  values  of  c  as  high  as  0.35.  In 
Figure  8  the  attenuation  constant  is  plotted  versus  ka 
for  a  fixed  concentration  c  =  0.209  using  the  SSA, 
WSA,  and  the  P-YA.  Values  for  the  WSA  for 
ka  <  0.75  are  not  shown  since  the  solution  fails  (Im 
K  <  0)  in  this  region.  However,  as  ka  increases,  it 
appears  that  the  WSA  tends  to  merge  with  the  P-YA. 
The  SSA,  on  the  other  hand,  consistently  overesti¬ 
mates  the  attenuation  over  the  full  range  of  ka 
values. 

4.  CONCLUSIONS 

This  paper  analyzes  the  effects  of  multiple  scatter¬ 
ing  on  the  coherent  wave  as  it  propagates  through  a 
discrete  random  medium  consisting  of  pair- 
correlated  scatterers.  At  low  concentrations 
(c  <  1%),  multiple  scattering  effects  are  seen  to  be 
negligible,  while  at  higher  concentrations,  suitable 
pair  correlation  models  must  be  assumed.  Compu¬ 
tations  are  presented  using  the  well-stirred  approxi¬ 
mation  and  the  Percus-Yevick  approximation  as  well 
as  a  model  due  to  Matern.  These  computations  are 
compared  (and  shown  to  be  in  excellent  agreement) 
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1  ig  1  The  coherent  attenuation  Im  l AC  ti  versus  concen¬ 
tration  i  at  t.i  -0  56  for  latex  spheres  in  water  using  viilTerent 
models  of  pair  correlation  functions 
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ka 

Fig.  8.  The  coherent  attenuation  Im  (K  k)  versus  ka  for  c  =  0.209 
for  late*  spheres  in  water  using  SSA,  P-YA.  and  WSA. 


with  recent  experiments  by  Ishimaru  (person  .  com¬ 
munication,  1981)  at  high  scatterer  concentrations 
(c  =  10%).  At  very  low  concentrations,  c  <,  1  %,  the 
WSA  is  compared  (and  shown  to  be  in  excellent 
agreement)  with  the  so-called  single  scattering  solu¬ 
tion  for  a  hypothetical  monodisperse  rain  medium 
(10* 4  <  c  <  10 ~2).  At  these  concentration  levels  the 
rain  medium  is  sparse  and  the  scatterers  are  essen¬ 
tially  decorrelated  so  that  multiple  scattering  effects, 
as  far  as  the  coherent  wave  attenuation  is  concerned, 
can  be  safely  neglected. 
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The  Effects  on  Pair  Correlation  Function  of  Coherent 
Wave  Attenuation  in  Discrete  Random  Media 

V.  N.  BRINGI,  MEMBER.  IEEE,  V.  V.  VARADAN,  AND 
V.  K.  VARADAN 

Abstract— The  Ptrcut-Yevick  approximation  (P-YA)  of  pair  cor¬ 
relation  function  for  hard  spheres  is  combined  with  the  7- matrix 
formulation  to  stud;  the  coherent  wave  attenuation  of  electromagnetic 
wave  propagation  in  a  discrete  random  medium.  The  effect  of  the  pair 
correlation  function  is  seen  to  be  significant  it  high  fractional  volumes 
of  the  discrete  scatterers  ( M.  12$),  but  at*  -  spends  on  the  frequency 
of  the  propagating  wave— the  effect  being  less  at  higher  frequencies. 
The  results  are  compared  with  previous  calculations  which  employed 
the  “well-stirred  approximation”  (WSA)  for  the  pair  correlation. 

INTRODUCTION 

We  consider  the  multiple  scattering  of  electromagnetic 
waves  by  randomly  distributed  dielectric  scatterers  using  the 
approach  given  by  Varadan,  Bringi,  and  Varadan  [  1  ] .  This 
formulation  uses  the  F-matrix  (2)  to  characterize  the  scat¬ 
tering  by  a  single  isolated  particle  followed  by  configurational 
averaging  techniques  (31,  [4|.  Lax’s  quasi-crystalline  approxi¬ 
mation  (QCA)  (4)  is  used  to  truncate  the  resulting  heirarchy 
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of  equations.  Thus,  only  knowledge  of  the  two  particle  (pair) 
correlation  function  is  required.  In  [  1 1 ,  we  assumed  that  the 
particles  were  hard  (nonpenetrating)  but  otherwise  uncor¬ 
related.  Talbot  and  Willis  [SI  refer  to  this  as  the  “well-stirred 
approximation”  (WSA).  This  yielded  a  set  of  “hole  correc¬ 
tion"  integrals  which  were  evaluated  analytically  and  the  ex¬ 
tinction  theorem  was  invoked  to  yield  the  dispersion  relations 
characterizing  the  bulk  or  effective  properties  of  the  medium 
which  was  solved  numerically  (1],  [6],  (7).  Computations 
of  the  effective  coherent  wave  attenuation  as  a  function  of 
frequency  were  presented  in  [1]  for  spherical  and  oblate 
spheroidal  scatterers  at  concentrations  (c)  of  0.05,  0.1,  and 
0.2.  The  WSA  leads  to  unphysical  nulls  in  the  plots  of  coher¬ 
ent  attenuation  as  a  function  of  frequency  or  concentration 
for  c  >  0. 1 25  and  in  fact  begins  to  fail  for  c  >  0.05  in  the  low 
frequency  or  Rayleigh  limit.  These  nulls,  however,  disappear 
at  higher  values  of  the  nondimensional  wave  number  ka  (a 
being  a  characteristic  dimension  of  the  obstacle). 

In  this  communication,  we  wish  to  consider  a  more  realis¬ 
tic  model  of  the  pair  correlation  function  which  is  valid  for 
higher  values  of  concentration.  The  Percus-Yevick  (P-YA) 
model  [8]  seems  most  suitable  at  the  present  time.  Calcula¬ 
tions  employing  P-YA  here  are  compared  with  WSA  results 
in  ( l  ] .  The  unphysical  nulls  in  [  1  ]  disappear  for  lower  values 
of  ka,  and  at  higher  values  of  ka  there  is  agreement  with  the 
results  obtained  in  [  1  ]  and  ( 7] . 

MULTIPLE  SCATTERING  FORMULATION 

We  consider  N(N  -*■  °°)  rotationally  symmetric  oriented 
scatterers  randomly  distributed  in  a  volume  V(V  -*  °°)  so  that 
the  number  of  particles  per  unit  volume  n0  =  N/V  is  finite, 
see  Fig.  1 .  Only  the  most  important  details  that  lead  to  the  dis¬ 
persion  equation  involving  the  pair  correlation  function  are 
presented  and  for  all  intermediate  steps,  we  refer  the  reader 
toll). 

Monochromatic,  plane  electromagnetic  waves  giving  rise  to 
an  electric  field  E°  are  assumed  to  propagate  parallel  to  the 
rotational  axis  of  symmetry  of  the  scatterers  (the  z-axis), 
see  Fig.  I .  The  field  scattered  by  the  ith  scatterer  is  denoted 
by  £j’  so  that  the  total  field  £  at  a  point  r  outside  the  scat¬ 
terers  is  given  by 

N 

E(?)  =  1°(?)  +  2  Et\n  C) 

i=l 

The  field  exciting  the  ith  scatterer  £*  is  given  by 
N 

£,'(?)  =  E°(?)  +  2  2/ V):  a<\r-r,\<2a.  (2) 

/*/ 

From  (1)  and  (2),  we  note  that 

£(?)  =  £,'(?)  +  *.'(?)  (3) 

so  that  the  exciting  and  scattered  fields  must  be  defined  in  a 
self-consistent  manner.  These  fields  are  expanded  in  a  set  of 
vector  spherical  functions  as  follows; 

•  /  o 

lV)='Z  2  2  Mm,  Re3WF-;,) 

1=0  m  =  0  o-t 

+  c‘omi  Re^om^-?.)};  a<\r-r,\<2a  (4) 
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Fig.  1.  Geometry  of  randomly  distributed  and  aligned  scattereis  with 
dielectric  constant  tr  excited  by  an  incident  plane  wave  £“(?).  Each 
scatterer  is  assumed  to  be  circumscribed  by  a  spherical  shell  of 
radius  a. 


Et\T)  =  2  2  2  {flo,m/^om/(r-?,) 

I  m  a 

+  cUlNon,tf-h)Y.  \r-rt\>2a  (5) 

where  {fc,  c}  and  (B,  C}  are  expansion  coefficients  of  the 
exciting  and  scattered  fields,  respectively.  The  vector  spherical 
functions  {A,  jV}  are  defined  by  Stratton  ( 9] . 

These  expansions  are  substituted  into  (2)  with  the  follow¬ 
ing  definition  of  the  T-matrix  of  a  single  scatterer 

r bL,  nr  (TWt)"  vx-tf*  i 

[c'm, 

•  P/  m  ']  (6) 

LVml'J 

(where  the  T-matrix  is  independent  of  the  position  of  the  scat- 
terers)  which  results  in  an  equation  for  the  exciting  field 
coefficients  {b,  c}  alone.  This  equation  is  averaged  over  the 
position  of  all  scatterers  where  the  QCA  [4)  is  involved, 
and  we  arrive  at  an  equation  for  the  configurational  average 
(A1) i  and  <c%  of  the  exciting  field  coefficients  with  one  par¬ 
ticle  held  fixed. 

We  assume  that  this  average  field  (the  coherent  field)  pro¬ 
pagates  in  a  medium  with  an  effective  complex  wavenumber 
a  —  (JC|  +  iK2)k  in  the  direction  of  the  original  incident  field 
in  the  discrete  medium.  Thus,  we  obtain 

(7) 


field.  The  dispersion  equations  then  take  the  following  form: 


=  2  2  2  «,(-«)xWx{roip[(r8!;)“ 

n=  1  p  =  1  X=ln-n'l 

•  *oo(n.  »\  X)  +  (rSip)2 ‘)Uo(». X)] 

+  ZeiPl(TVi^Mn,n',  A) 

+  (rM?)22x*o("."'.x)]} 

in'Ztln- 

o*  M  n+n' 

-SS  2.  i^-i)Vtf)x{l'oip[(7'ol2>,‘ 

n=l  p=l  X=|n-n’l 

*  Xo eO*.  n'<  X)  +  (rSlp)2*  X)] 

+  2elpt(r?«;),2Xoe(«,*»'.X) 

+  (n{?)22^e(».»'.x)]}  ( 

where  (///)*  is  an  integral  given  by 


{JH)x(K,k,c) 


(ka)2-(Ka? 


[2kajx(2Ka)hK\2ka) 


(cUlh=i,Zomlel*'*< 


-  2Kahx(2ka)jx\2Ka)] 

+  24c  f  x2  [*(*)  -  1  ]  hx(kx)jK(Kx)  dx 


and  the  functions  {i^,  x)  uc  defined  in  (7). 

In  the  above  equation,  c  -  4ira3n0/3  is  the  effective  “spher¬ 
ical”  concentration.  In  (9)  and  (10),  \p  and  x  are  independent 
of  k  and  K,  and  expressions  for  them  may  be  found  in  [  1 )  and 
(7).  In  (11),  and  hx  are  the  spherical  Bessel  and  Hankel 
functions  and  g(x)  is  the  pair  correlation  function  that  is  dis¬ 
cussed  below. 

PERCUS-YEVICK  APPROXIMATION  FOR  THE 
PAIR  CORRELATION  FUNCTION 

It  is  well  known  that  at  high  concentrations  the  effects,  of 
pair  correlation  become  important.  Following  the  notation  in 
(1),  the  conditional  probability  distribution  p(fj  | ?,)  is  ex¬ 
pressed  as 


„  l—  g(x);  x  >  l 
P(?i  I  '<)  =  ]  V 


where  {  Y,  Z }  are  expansion  coefficients  of  the  average  exciting  where  x  -  ]ft  —  rf  1/2 a.  The  Percus-Yevick  integral  equation 
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[8]  for  hard  spheres  is  given  by 


t(2)=1+*0  /  T(x')dx' 

—  n0  /  T(x')r(2  —  x')dx'  (13) 

I  3 '  l<2« 

I  i?  —  J'l>2a 


The  function  r(x)  is  related  to  the  direct  correlation  function 
C(x)  of  Ornstein  and  Zemike  and  the  pair  distribution  func¬ 
tion  g(x)[9): 


g(x)  —  f(x); 

|x|>  2a 

*(i)  ~  0; 

|x  |  <  2j 

C(J)  =  r(x); 

|x|  <  2j 

C(x)  =  0; 

|i|>2u. 

(14) 


For  hard  spheres,  the  direct  correlation  function  C(x )  is 
known  explicitly  and  Wertheim  [10)  has  obtained  a  series 
solution  for  g(x)  as  a  function  of  concentration  when  g(x) 
is  radially  symmetric.  Throop  and  Bearman  [11]  have  used 
the  Wertheim  result  and  provided  tabulated  values  of  g(x)  as  a 
function  of  x  for  several  values  of  c. 

The  highest  concentration  for  which  the  tabulated  results 
can  be  used  is  c  =  0.26.  Beyond  this  concentration  g(x) 
oscillates  significantly  from  its  asymptotic  value  of  1  for  x  > 
4.  Thus,  the  integral  in  (1 1)  cannot  be  evaluated  accurately. 

At  low  values  c,  g(x)  =*  1  and  hence  the  integral  in  (1 1)  is 
negligible  and  the  remaining  term  is  simply  the  WSA  which 
was  used  in  [l].  Thus,  (11)  can  be  regarded  as  a  modified 
“hole  correction  integral'’  and  is  of  the  same  form  used  by 
Twersky  [  12] ,  [  13) .  Talbot  and  Willis  [5]  have  also  suggested 
models  of  the  pair  correlation  function  given  by  Matem  which 
are  valid  for  c  <  0.125.  The  advantage  of  the  Matem  model  is 
that  it  is  completely  analytic. 


Fig.  2.  Coherent  wave  attenuation,  Im  (K/k)  as  a  function  of  concen¬ 
tration  for  dielectric  spheres. 


Fig.  3.  Attenuation  coefficient  a  as  a  function  of  ka  for  spherical  di¬ 
electric  scatterers. 


NUMERICAL  COMPUTATIONS 

In  the  Rayleigh  limit,  Twersky  (13)  has  given  an  expression 
for  Ki/k  by  considering  the  leading  effects  of  the  pair  correla¬ 
tion  function: 


K2/k  =  Im  (K/k)  =  c(ka)3 


er-  1 

er+2 


W 


(15) 


where  e,  is  the  dielectric  constant  of  the  scatterer  and  W  is  the 
packing  factor  given  by 

(i-c)4  r  , 

W  = - —  =  1  +  24c  /  x2[g(x)- lj  dx.  (16) 

(1  +  2c)2  /0 

In  Fig.  2,  we  compare  Twersky’s  result  for  the  coherent 
attenuation  with  computations  for  ka  =  0.05  and  er  =  3.168 
as  a  function  of  concentration.  Also  shown  is  the  result  assum¬ 
ing  the  WSA.  It  is  clear  that  the  WSA  fails  for  c  >  0.05  while 
the  present  calculations  is  in  good  agreement  with  (15). 

In  Figs.  3,  4,  and  5,  calculations  using  the  P-YA  model  are 
presented  as  a  function  of  ka  for  c  =  0.2  for  spheres  and 


Fig.  4.  Attenuation  coefficient  a  versus  ka  for  oblate  spheroidal 
scatterers. 
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Fig.  5.  Attenuation  coefficient  a  versus  ko  for  oblate  spheroidal 
scatterers. 

oblate  spheroids.  Equations  (9)  and  (10)  were  solved  numeri¬ 
cally  to  find  K(k)  using  analytical  results  at  low  frequencies 
(Rayleigh  limit)  as  initial  guesses  in  a  root  searching  algorithm. 
For  purposes  of  comparison,  [1.  Figs.  3,  3,  and  6]  using  just 
the  WSA,  are  reproduced.  It  is  clear  that  the  nulls  appearing 
as  a  result  of  the  WSA  are  wrong  and  unphysicaL  However, 
it  is  interesting  to  observe  that  for  ka  >  1.5,  it  appears  that  the 
results  obtained  using  the  WSA  alone  tend  to  the  P-YA  result. 
Note  that  in  ( 1,  Fig.  3]  the  curves  marked  c  =  0.05  and  c  = 
0. 1  must  be  interchanged. 

CONCLUSION 
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Calculations  using  the  Percus-Yevick  (P-YA)  pair  coiTela- 
•V*  tion  function  for  hard  spheres  are  compared  with  previous 
computations  [1]  using  the  well  stirred  approximation.  The 
effects  of  pair  correlation  are  seen  to  be  significant  for  high 
|  values  of  scatterer  concentrations  c  >  0.05,  although  this 
effect  seems  to  decrease  with  increasing  frequency.  It  is  clear 
that  the  WSA  cannot  be  used  for  arbitrary  concentrations  and 
some  of  the  results  presented  in  [1],  [6],  and  [14]-[16]  are 
thus  unphysical. 
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Attenuation  of  electromagnetic  waves  by  a  random  distribution  of  pair-correlated  dielectric  spheres 
is  studied  as  a  function  of  frequency  and  volume  concentration  of  spheres.  The  main  aim  of  this  paper 
is  to  compare  theoretical  results  obtained  using  a  self-consistent  multiple-scattering  formulation  and 
measured  values  of  attenuation  for  latex  spheres  in  water.  The  agreement  between  theory  and  experi¬ 
ment  is  very  good. 


INTRODUCTION 

A  discrete  random  medium  is  defined  here  as  a 
random  distribution  of  a  large  number  of  identical 
pair-correlated  scatterers  embedded  in  an  infinite  ho¬ 
mogeneous  matrix  medium.  The  fractional  volume 
concentration  c  is  assumed  uniform,  and  the  effects  of 
pair  correlation  are  described  by  the  radial  distri¬ 
bution  functions  arising  in  the  statistical-mechanical 
treatment  of  dense  gases  and  liquids.  We  consider 
plane  wave  propagation  through  such  a  model 
random  medium  with  attention  being  focused  on  the 
coherent  (or  ensemble  averaged)  wave  attenuation. 
The  discrete  random  medium  is  then  characterized 
by  an  effective,  complex  wave  number  K ,  or  equiva¬ 
lently,  by  an  effective  dielectric  constant.  If  the  scat¬ 
terers  and  the  matrix  medium  are  assumed  lossless, 
then  the  coherent  wave  attenuation  is  caused  by  scat¬ 
tering,  and  it  is  of  great  interest  to  determine  the 
wave  frequency  and  scatterer  concentration  for 
which  multiple-scattering  effects  dominate.  Extensive 
work  by  Twersky  [1977,  1978a,  b,  c,  1982]  has  laid 
the  foundation  for  multiple-scattering  theory  in  dis¬ 
crete  random  media.  A  related  approach  using  the  T 
matrix  of  a  single  scatterer  [ Varadan  and  Varadan, 
1980a]  together  with  configurational  averaging  pro¬ 

1  Wave  Propagation  Group,  Boyd  Laboratory,  Department 
of  Engineering  Mechanics,  and  Atmospheric  Sciences  Program, 
Ohio  State  University,  Columbus,  Ohio  43210. 

2  Department  of  Electrical  Engineering,  Colorado  State  Uni¬ 
versity,  Ft.  Collins,  Colorado  80523. 

*  Department  of  Electrical  Engineering,  University  of  Wash¬ 
ington.  Seattle.  Washington  98195. 

Copyright  1983  by  the  American  Geophysical  Union. 

Paper  number  3S0089. 

0048-6604  83  0506-0089S08.00 


cedures  has  been  used  by  the  authors  to  develop  a 
computational  method  for  handling  acoustic,  electro¬ 
magnetic,  and  elastic  waves  [ Varadan  et  a/.,  1979, 
1982;  Varadan  and  Varadan,  1980a,  b;  Bringi  et  a/„ 
1981,  1982a,  b],  Lax’s  [1952]  quasi-crystalline  ap¬ 
proximation  (QCA)  is  used  in  conjunction  with 
various  models  for  the  radial  distribution  function 
enabling  passage  to  an  average  dielectric  medium 
whose  properties  depend  on  frequency,  scatterer  con¬ 
centration,  and  the  dielectric  constants  of  both  the 
scatterers  and  the  matrix  medium.  The  intent  of  this 
paper  is  to  compare  theoretical  calculations  of  coher¬ 
ent  attenuation  with  optical  experiments  conducted 
recently  by  Ishimaru  and  Kuga  [1982],  Preliminary 
comparisons  between  theory  and  experiment  have 
been  reported  [ Bringi  et  al.,  1982 a,  b],  where  the 
single-scattering  approximation  is  compared  with 
multiple-scattering  calculations  using  various  forms 
for  the  radial  distribution  function  (e.g.,  Matern  and 
Percus-Yevick  models). 

In  this  paper  we  provide  an  outline  of  the 
multiple-scattering  theory  followed  by  a  brief  dis¬ 
cussion  of  the  self-consistent  radial  distribution  func¬ 
tion  used  in  the  computations.  This  is  followed  by 
calculation  of  coherent  wave  attenuation  as  a  func¬ 
tion  of  scatterer  concentration  (0  <  c  <  0.4)  for  a 
number  of  wave  frequencies  (0  <  ka  <  7)  coinciding 
with  the  experimental  setup  of  Ishimaru  and  Kuga 
[1982],  Finally,  previous  computations  presented  by 
Bringi  et  al.  [1981]  corresponding  to  the  experiments 
of  Hawley  et  al.  [1967]  are  discussed  in  the  light  of 
the  results  presented  here. 

OUTLINE  OF  THEORY 

Consider  an  incident  TEM  wave  propagating  in  an 
infinite  lossless,  background  medium  of  cm,  g0  which 
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is  filled  with  a  random  distribution  of  identical 
dielectric  spheres  of  e,  n0 .  The  total  dielectric  field  at 
any  point  in  the  background  medium  is  the  sum  of 
the  incident  field  and  the  fields  scattered  by  all  the 
scatterers, 

H 

E(r)  =  E‘«(r)+  XEJfr-rJ  (I) 

i=  I 

where  E*(r  —  rj  is  the  field  scattered  by  the  ith  scat- 
terer  at  the  observation  point  r  and  r,  is  the  position 
of  the  ith  scatterer.  The  field  that  excites  the  ith  scat¬ 
terer,  however,  is  the  incident  field  plus  the  fields 
scattered  from  all  the  other  scatterers.  The  term  ex¬ 
citing  field,  E*.  is  used  to  distinguish  between  the  field 
actually  incident  on  a  scatterer  and  the  external  inci¬ 
dent  field  Ela<  produced  by  a  source  at  infinity.  Thus 
at  a  point  r  in  the  vicinity  of  the  ith  scatterer,  the 
exciting  field  is  expressed  as 

,V 

E*(r)  =  E'"c(r)  +  £  E*(r  -  tj)  (2) 

i*  i 

From  (1)  and  (2)  we  note  that 

E(r)  -  E*(r)  +  EJ(r)  (3) 

so  that  the  exciting  and  scattered  fields  must  be  de¬ 
fined  in  a  self-consistent  manner.  Waterman’s  T 
matrix  formalism  is  used  to  relate  the  exciting  and 
scattered  field  expansion  coefficients  (the  fields  are 
expanded  using  the  vector  spherical  functions  as  a 
basis).  Configurational  averaging  is  performed  in  (2) 
as  described  by  Varadan  et  al.  [1979].  The  T  matrix 
together  with  the  QCA  is  used  to  generate  a  homoge¬ 
neous  system  of  equations  whose  singular  solutions 
yield  the  average  propagation  constant  (K  =  + 

jK2)  for  the  discrete  random  medium.  The  system  of 
equations  is  given  below,  and  the  reader  is  referred  to 
Bringi  et  al.  [1981],  Varadan  et  al.  [1979,  1982],  and 
Varadan  and  Varadan  [1980a]  for  full  details. 

XX  X 

rm  =  Z  I  I  (JH), 

«■  1  1  |»-tw| 

[Y,  TljlMn,  m,  ;.)  +  Z,  m,  /.)]  (4a) 

XX  X 

Z.  -  l  I  I  IWi 

■■  i  f  i  i* i«-«i 

[-  Y,T"x(ii,  m,  a)  +  Z'TjtjMn,  m,  /.)]  (46) 

where 

{JH)k  “  (ko)1  -  (Ka)1 

[2ka/,,(2Ka)6i(2ka)  —  2Kah^2ka\j'i{2Ka)'\ 


+  24c  J  xJ{g(x)  -  l}6i(2/cax)/i(26Cax)  dx  (4c) 

In  (4a)  and  (46),  Ym  and  Zm  are  the  unknown  ampli¬ 
tudes  of  the  average  exciting  field,  T„p  refers  to  the 
elements  of  the  T  matrix  for  a  sphere,  and  the  func¬ 
tions  <1/  and  x  are  defined  by  Bringi  et  al.  [1981].  In 
(4c),  g(x)  is  the  radial  distribution  function;  jx  and  hx 
are  the  spherical  Bessel  functions,  and  the  primes 
denote  differentiation  with  respect  to  the  argument;  k 
is  the  wave  number  of  the  background  medium;  a  is 
the  scatterer  radius;  and  c  is  the  fractional  con¬ 
centration  of  the  scatterers.  Assuming  that  ^(x),  ka ,  c, 
and  the  T  matrix  are  known,  the  singular  value  of 
the  coefficient  matrix  generated  from  (4)  can  be 
solved  for  the  average  propagation  constant  K  = 
K  i  +jKt. 

RADIAL  DISTRIBUTION  FUNCTION 

The  discrete  random  medium  is  considered  as  a 
statistical  ensemble  of  impenetrable  spheres.  In  the 
statistical  mechanics  literature  this  is  synonymous 
with  an  ensemble  of ‘hard’  spheres.  The  radial  distri¬ 
bution  g(r)  is  defined  in  terms  of  the  two-particle 
joint  probability  density  p( r2 1  r()  defined  as 

Pbjltd^-gdrj-rM  |ry  — r(|s2  a 

(5) 

Pfrj|r,)  =  0  I r>  —  r( |  <  2a 

Equation  (5)  implies  that  the  particles  are  hard  (no 
interpenetration)  and  the  excluded  volume  is  a  sphere 
of  radius  a  although  the  particles  themselves  may  be 
nonspherical.  The  function  g{r),  (r  =  r0),  is  called  the 
pair  correlation  function  and  depends  only  on 
i Tij I  -  |ry  —  r,|  because  of  translational  invariance  of 
the  system  under  consideration.  Several  theories  and 
calculations  are  available  for  determining  g(r), 
namely,  the  hypernetted-chain  equation,  the  Percus- 
Yevick  approximation  (P-YA),  the  self-consistent  ap¬ 
proximation,  Monte  Carlo  calculations,  etc.  These 
approaches  are  described  by  McQuarrie  [1976],  but 
the  reader  is  referred  to  any  standard  text  on  statis¬ 
tical  mechanics  for  full  details. 

The  pair  correlation  function  for  an  ensemble  of 
particles  depends  on  the  nature  and  range  of  the  in¬ 
terparticle  forces.  The  average  of  several  measure¬ 
ments  of  a  statistical  variable  that  characterizes  an 
ensemble  will  depend  on  the  pair  correlation  func¬ 
tion.  To  obtain  expressions  for  the  pair  correlation 
function,  one  needs  a  description  of  the  interparticle 
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forces.  In  our  case,  we  assume  that  the  scatterers 
behave  like  effective  hard  (impenetrable)  spheres. 
Percus  and  Yevick  [1958]  have  obtained  an  approxi¬ 
mate  integral  equation  for  the  pair  correlation  func¬ 
tion  of  a  classical  fluid  in  equilibrium.  Wertheim 
[1963]  has  obtained  a  series  solution  of  the  integral 
equation  for  an  ensemble  of  hard  spheres.  The  statis¬ 
tics  of  the  fluid  are  then  the  same  as  those  of  the 
ensemble  of  discrete  hard  particles  that  we  are  con¬ 
sidering. 

Although  integral  expressions  for  the  correlation 
functions  also  result  in  a  hierarchy,  Percus  and  Yevick 
[1958]  have  truncated  the  hierarchy  by  making  cer¬ 
tain  approximations  that  result  in  a  self-consistent 
relation  between  the  pair  correlation  function  g(r) 
and  the  direct  correlation  function  C(r).  The  direct 
correlation  can  be  defined  to  be  the  direct  effect  of 
scatterer  1  on  2  which  roughly  has  the  range  of  the 
interparticle  potential  and  also  an  indirect  effect  due 
to  the  effect  of  1  on  3  which  in  turn  has  an  effect  on 
2.  Since  3  can  be  any  scatterer  of  the  system,  it  will 
include  a  sum  on  all  scatterers  and  an  integration 
over  the  volume  of  the  system.  Fisher  [1965]  com¬ 
ments  that  the  P-YA  is  a  strong  statement  of  the 
extremely  short  range  nature  of  the  direct  correlation 
function.  For  impenetrable  spheres,  the  range  of  the 
interparticle  potential  is  2a.  The  integral  equation 
has  the  form  [see  Percus  and  Yevick,  1958]  given  by 

r(r)  =  1  +  «0  r(r  )  dr'  -  n0  r(r')r(r  -  r)  dr' 

Jr<2a  J 

r><2a 
|r  -r»|  >  la 

(6) 

where 


r(r)  =  c/(r) 

r  >  2a 

gir)  =  0 

r  <  2a 

r(r)=  —  C(r) 

r  <  2a 

C{r)  =  0 

r  >  2a 

Wertheim  [1963]  has  solved  the  integral  equation  by 
Laplace  transformation  that  results  in  an  analytic  ex¬ 
pression  for  C(r)  in  the  form 

C(r)  =  -(1  -  »j)  "‘[(1  +  2g)2  -  6>jU  +  ]»j)2r 

+  *1  +  2*))V  2]  >j  =  c  8  (8) 

where  c  is  the  effective  spherical  concentration  of  the 
particles.  The  P-YA  fails  as  the  concentration  ap¬ 
proaches  the  close  packing  factor  for  spheres  and  is 
expected  to  be  good  for  c  <  0.3  or  0.4. 


Fig.  1.  The  Percus-Yevick  (P-YA)  pair  correlation  function  for 
‘hard  ‘  spheres. 


Equation  (8)  can  be  substituted  back  into  (6)  to 
yield  a  series  solution  for  g(r)  in  the  form  [see 
Wertheim,  1963] 

3C 

£/(r)=  £  0„(r)  (9) 

n  -  l 

where 

g.(r)  =  ^  |  e',r  "’[L(f)|S(t)]"r  dt  (10) 

S(t)  =  (I  -  tj2)rJ  +  6>;(1  -  t))r2  +  18t/2f  -  12^(1  +  2r|)  (II) 

Ut)  =  12tj[(l  +r,2)t  +(1  +2*1]  (12) 

Throop  and  Bearman  [1965]  have  tabulated  t/(r)  as  a 
function  of  r  for  values  of  r\  =  c/8.  A  few  repre¬ 
sentative  plots  of  the  Percus-Yevick  (P-Y)  pair  corre¬ 
lation  function  are  shown  in  Figure  1. 

Another  approximation  to  cj(r)  is  the  hypernetted- 
chain  equation  (HNC)  which  differs  from  the  P-YA 
in  that  Ihe  direct  correlation  function  C(r)  has  a 
longer  range.  In  general,  P-YA  is  expected  to  be 
better  than  HNC.  The  P-Y  and  HNC  results  are 
good  approximations  to  </(r)  at  low  concentrations 
but  are  appreciably  in  error  at  high  concentrations. 

At  low  concentrations,  a  series  expansion  can  be 
obtained  from  (6)  by  iterating  in  powers  of  n0.  This 
power  series  in  density  is  the  virial  expansion  and  has 
been  used  at  low  frequencies  by  Twersky  [1977, 
1978a,  h,  c-],  Bringi  et  at.  [1981,  1982a.  h],  and 
Varadan  et  al.  [1982],  For  the  leading  terms,  (6)  is 
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Fig.  2.  The  P-YA,  hole  correction,  virial  expansion  pair  corre¬ 
lation  functions,  and  pair  correlation  functions  based  on  module 
histogram  for  impenetrable  spheres. 


written  as  [see  Percus  and  Yevick ,  1958] 
*12  =  1  +  "o  J/2J  T23(T,3  -  I)  dr3 

+  n0  J* /23/13  *23  rl3  dr 3 


fj~~  1  I  r,>  |  <  2a 


/«-  0 


|r,yl>  2fl 


Commencing  at  rl2  =  1,  one  can  readily  iterate  a 
power  series  to  obtain 


rw  =  I  + 


”0  J* /23/13  dr 3  +  n*  j*j* [/2 


+  ~fn  fa  dtl  </r4  (14) 

For  hard  spheres,  g(r)  can  be  determined  exactly  to 
0(n 0)  [see  Percus  and  Yevick,  1958;  Twersky,  1977, 
1978a,  b.  c]  given  by 


y{r)  =1  r  >  2a 


4tt  ,  (  3r  r3  \ 


<  r  <  2a  (15) 


g{r)  =  0  r  <  a 


but  a  slight  discrepancy  has  been  observed  in  0(n£) 
[see  McQuarrie,  1976]  because  of  the  approximation 
in  the  P-YA.  The  variation  of  g(r)  in  terms  of  virial 
expansion  is  shown  in  Figure  2.  At  low  frequencies,  a 
relatively  simpler  form  of  g(r)  given  in  terms  of  trig¬ 
onometrical  functions  as  shown  in  Figure  2  may  also 
be  used. 

At  higher  concentrations,  the  values  of  </(r)  using 


the  P-YA  are  in  appreciable  error.  The  reason  for 
this  error  is  that  when  an  approximate  theory  such 
as  P-YA  for  g(r )  is  used  to  calculate  the  pressure,  two 
different  values  are  obtained  from  the  pressure  equa¬ 
tion  and  the  compressibility  equation; 


=  K~T  {”°  ~  | V  ,  K‘T  ~  *  *  "  T(r'  ”  *  df 


P  =  n0K 


r(r)r  •  ye~V{,)K,T  i 


where  KB  is  the  Boltzmann  constant  and  T  is  the 
system  temperature.  For  hard  spheres. 


K(r)  =oo  |r|  <  2a 
F(r)  =  0  |  r  |  >  2a 


(17) 


which  shows  that  the  system  is  independent  of  tem¬ 
perature  but  depends  only  on  number  density  or  con¬ 
centration.  An  exact  theory  for  g(r)  should  give  the 
same  calculated  pressure  from  both  equations.  Row- 
finson  [1965]  has  suggested  a  method  known  as  the 
self-consistent  approximation  (SCA)  for  optimizing 
the  P-YA  and  removing  the  inconsistency  in  the  two 
pressure  equations  (16)  by  assuming  that  the  direct 
correlation  function  may  be  written  as 


^scaW  -  Cp.vAlr)  +  d>CHNC(r)  (18) 


where  d>  is  an  adjustable  parameter  which  depends 
on  concentration  but  not  on  separation  r  and  is 
the  self-consistent  approximation  of  the  direct  corre¬ 
lation  function.  Using  (7)  and  (18),  the  two  pressure 
equations  (16)  that  depend  on  g(r )  are  solved  self- 
consistently  by  adjusting  the  function  <t>  and  lead  to 
an  integral  equation  for  g{r)  which  is  then  solved 
numerically.  At  higher  concentrations,  the  SCA  is 
thus  an  improvement  over  the  P-YA  and  HNC  ap¬ 
proximation.  From  (18)  it  is  obvious  that  when 
<t>  =  0,  (18)  reduces  to  the  P-YA.  whereas  if  d>  =  1,  it 
gives  the  HNC  approximation.  Reed  and  Gubbins 
[1973]  have  provided  tabulated  values  of  </scA(r)  for 
0.0524  ^  c  <,  0.417  and  for  1.0  <,  r  <;  5.0,  and  a  gen¬ 
eral  computer  program  for  computing  g^^  is  given 
by  McQuarrie  [1976].  Tabulated  values  of  </SCA(r)  are 
also  given  by  Ikdtrs  and  Henderson  [1969],  but  they 
are  not  useful  for  computations  because  the  values  of 
g  are  given  only  for  a  short  range  of  r.  Monte  Carlo 
calculations  for  g(r )  are  done,  and  their  values  are 
tabulated  by  Barker  and  Henderson  [1971];  they 
cannot  be  used  for  computations  for  the  same  rea¬ 
sons.  However,  (he  values  given  by  Barker  and  Hen¬ 
derson  [1971]  provide  a  very  good  check  for  higher- 
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Fig.  3.  The  self-consistent  pair  correlation  function  for  ‘hard’ 
spheres. 

order  virial  coefficients.  The  computer  program  given 
in  the  work  by  McQuarrie  [1976]  is  thus  preferable 
for  computational  purposes.  A  representative  plot  of 
0scA(r)  versus  r  is  shown  in  Figure  3. 

COMPARISON  WITH  EXPERIMENTS 

Extensive  propagation  experiments  have  been  con¬ 
ducted  by  Ishimaru  and  Kuga  [1982]  where  the  dis¬ 
crete  random  medium  comprises  a  collection  of 
nearly  identical  latex  spheres  embedded  in  water  and 
illuminated  with  a  HeNe  laser  (/  =  0.6328  m).  The 
coherent  wave  attenuation  is  measured  over  a  range 
of  concentrations  (up  to  40%)  and  for  a  number  of 
equivalent  particle  sizes  (low.  Mie,  and  optical  ka 


Fig.  4  Im  (K  It)  versus  concentration  using  the  self-consistent 
radial  distribution  function. 


Fig.  5.  Im  (K,k)  versus  concentration  using  the  self-consistent 
radial  distribution  function. 


values).  Computations  of  coherent  attenuation  K/k 
as  a  function  of  frequency  (or  ka)  and  concentration 
are  performed  by  determining  the  singular  values  of 
the  coefficient  matrix  obtained  from  (4a)  and  (46) 
using  the  self-consistent  radial  distribution  function. 
(K/k)1  is  the  effective  dielectric  constant  of  the  bulk 
random  medium  relative  to  the  background  water 
medium.  Figure  4  shows  computations  of  Im  (K/k)  as 
a  function  of  concentration  for  ka  values  of  0.529  and 
0.681  together  with  the  experimental  points  of  Ishi¬ 
maru  and  Kuga  [1982],  Note  that  the  original  experi¬ 
mental  points  have  been  converted  to  Im  (K/k)  by 
using  a  normalization  value  described  by  Ishimaru 
and  Kuga  [1982].  The  comparison  between  theory 
and  experiment  is  seen  to  be  excellent  even  at  high 
concentrations.  In  Figure  5  we  show  computations  of 
Im  (K/k)  versus  concentration  of  ka  values  of  3.518 
and  7.28.  Again  the  experimental  points  are  seen  to 
be  in  excellent  agreement  with  the  computations. 
Even  though  the  experiments  were  not  conducted  at 
c  >  10%  at  these  ka  values,  the  behavior  of  Im  (K'k) 
is  correctly  reproduced  by  the  self-consistent  radial 
distribution  function.  In  Figure  6  we  compare  our 
computations  with  the  experimental  results  of  Ishi¬ 
maru  and  Kuga  [1982].  The  parameter  y  =  2K3fn0o, 
where  K2  =  Im  (K).  n0  is  the  number  density,  and  <r, 
is  the  extinction  cross  section  of  a  single  sphere.  We 
remark  that  n0a,  is  the  approximate  value  of  the 
attenuation  at  very  low  concentrations.  Hence  the 
parameter  y  is  the  attenuation  at  any  concentration 
normalized  to  the  low  concentration  limit.  Thus,  as 
c-»0,  y  ♦  1.  Further,  as  explained  by  Ishimaru  and 
Kuga  [1982],  since  the  latex  spheres  used  in  the  ex- 
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Fig.  S.  Im  (K;k)  versus  concentration  using  the  self-consistent 
Fig.  3.  The  self-consistent  pair  correlation  function  for  'hard'  radial  distribution  function, 

spheres. 


order  virial  coefficients.  The  computer  program  given 
in  the  work  by  McQuarrie  [1976]  is  thus  preferable 
for  computational  purposes.  A  representative  plot  of 
<fecA(r)  versus  r  is  shown  in  Figure  3. 

COMPARISON  WITH  EXPERIMENTS 

Extensive  propagation  experiments  have  been  con¬ 
ducted  by  Ishimaru  and  Kuga  [1982]  where  the  dis¬ 
crete  random  medium  comprises  a  collection  of 
nearly  identical  latex  spheres  embedded  in  water  and 
illuminated  with  a  HeNe  laser  (/.  =  0.6328  m).  The 
coherent  wave  attenuation  is  measured  over  a  range 
of  concentrations  (up  to  40%)  and  for  a  number  of 
equivalent  particle  sizes  (low.  Mie,  and  optical  ka 


Fig  4.  Im  (fC  k)  versus  concentration  using  the  self-consistent 
radial  distribution  function. 


values).  Computations  of  coherent  attenuation  K/k 
as  a  function  of  frequency  (or  lea)  and  concentration 
are  performed  by  determining  the  singular  values  of 
the  coefficient  matrix  obtained  from  (4a)  and  (46) 
using  the  self-consistent  radial  distribution  function. 
(K/k)2  is  the  effective  dielectric  constant  of  the  bulk 
random  medium  relative  to  the  background  water 
medium.  Figure  4  shows  computations  of  Im  (K/k)  as 
a  function  of  concentration  for  ka  values  of  0.529  and 
0.681  together  with  the  experimental  points  of  Ishi¬ 
maru  and  Kuga  [1982].  Note  that  the  original  experi¬ 
mental  points  have  been  converted  to  Im  (K/k)  by 
using  a  normalization  value  described  by  Ishimaru 
and  Kuga  [1982].  The  comparison  between  theory 
and  experiment  is  seen  to  be  excellent  even  at  high 
concentrations.  In  Figure  5  we  show  computations  of 
Im  (K/k)  versus  concentration  of  ka  values  of  3.518 
and  7.28.  Again  the  experimental  points  are  seen  to 
be  in  excellent  agreement  with  the  computations. 
Even  though  the  experiments  were  not  conducted  at 
c  >  10%  at  these  ka  values,  the  behavior  of  Im  (AC  A.) 
is  correctly  reproduced  by  the  self-consistent  radial 
distribution  function.  In  Figure  6  we  compare  our 
computations  with  the  experimental  results  of  Ishi¬ 
maru  and  Kuga  [1982].  The  parameter  y  =  2K:  'n0o, 
where  AC2  =  Im  (AC),  n0  is  the  number  density,  and  <r, 
is  the  extinction  cross  section  of  a  single  sphere.  We 
remark  that  n0a,  is  the  approximate  value  of  the 
attenuation  at  very  low  concentrations.  Hence  the 
parameter  y  is  the  attenuation  at  any  concentration 
normalized  to  the  low  concentration  limit.  Thus,  as 
c-»0,  */■■■♦  1.  Further,  as  explained  by  Ishimaru  and 
Kuga  [1982],  since  the  latex  spheres  used  in  the  ex- 
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Fig.  6.  Plot  ofy  =  2  Im  K/n0a ,  versus  concentration  where  n0 
is  the  number  density  and  a,  is  the  extinction  cross  section  of  a 
single  sphere  with  relevant  ka. 


periment  were  not  exactly  identical,  the  measured 
value  of  the  attenuation  at  c  <,  0.0038  did  not  exactly 
correspond  to  n0a,  where  a,  was  computed  at  the 
average  measured  radius  of  the  sphere.  In  order  to 
take  care  of  this  discrepancy,  the  measured  value  of 
the  attenuation  i  at  low  values  of  c  was  used  to 
define  the  effective  radius  of  a  sphere  which  pre¬ 
dicts  the  correct  attenuation  in  the  single-scattering 
approximation.  All  experimental  results  were  nor¬ 
malized  with  respect  to  x0  =  n0  <r,(ntff)  and  y  =  x/x0 . 
Thus  in  all  comparisons  we  have  also  normalized  the 
computed  value  of  x  by  x0  using  the  effective  radius 
supplied  in  Table  I  of  Ishimaru  and  Kuga  [1982].  The 
computations  are  in  excellent  agreement  with  experi¬ 
ment  for  ka  values  of  3.518  and  7.28  (see  Figure  6). 
Note  the  sharp  decrease  in  y  for  c  >  10%,  especially 
for  the  ka  value  of  7.28.  The  dramatic  importance  of 
multiple  scattering  is  seen  clearly  from  Figure  6  for 
c  >  1%:  this  value  of  c.  however,  depends  on  ka. 

In  our  previous  paper  [Bringi  et  al.,  1982 a,  b ]  we 
observed  that  at  higher  values  of  ka,  a  less  rigor¬ 
ous  form  for  the  radial  distribution  function 
</(/•)  =  1,(1  -  c)  may  suffice  even  at  relatively  high 
concentrations  (c  >  40%).  It  is  important  to  note 
that  at  wavelengths  comparable  to  obstacle  size  and 
higher,  the  scattering  is  mostly  in  the  forward  direc¬ 
tion.  Thus,  in  this  case,  repeated  scattering  should 
not  be  important,  since  the  backscattered  wave  is 
significantly  smaller  than  the  forward  scattered  wave. 
This  may  explain  why  the  computations  shown  in 
Figure  7  for  ka  =  11.8  are  in  reasonable  agreement 
with  the  experiment  done  by  Hawley  et  al.  at  higher 
values  of  concentrations.  It  should  be  noted  that  the 
P-YA  and  SCA  fail  as  the  concentration  approaches 


Fig.  7.  Attenuation  constant  2  Im  {K'k)  versus  concentration 
for  Aa  =  1 1 .8  using  the  ‘well-stirred’  (effective  concentration) 
radial  distribution  function.  The  experimental  points  are  taken 
from  Hawley  et  al.  [1967], 

the  close  packing  factor  for  spheres  and  are  expected 
to  be  good  for  c  <  0.42%.  Thus,  based  on  the  results 
presented  in  Figure  7,  the  coherent  wave  attenuation 
is  less  sensitive  to  the  precise  form  of  the  pair  corre¬ 
lation  function  at  high  concentrations  for  relatively 
high  values  of  ka. 

CONCLUSIONS 

This  paper  has  highlighted  the  comparison  be¬ 
tween  theory  and  experiments  for  coherent  wave  at¬ 
tenuation  in  a  discrete  random  medium.  Multiple¬ 
scattering  effects  must  be  considered  for  volume  con¬ 
centrations  exceeding  around  1%  especially  for  ka 
values  of  <10.  At  higher  concentrations  (c  <1  40%) 
and  for  ka  <  10  the  form  of  the  radial  distribution 
function  is  very  important.  At  values  of  ka  >  10,  cal¬ 
culations  show  that  the  coherent  wave  attenuation  is 
less  sensitive  to  the  precise  form  of  the  radial  distri¬ 
bution  function.  Computations  and  experimental  re¬ 
sults  suggest  that  at  high  ka  values  the  simple  ‘well- 
stirred’  approximation  for  the  radial  distribution 
function  (</(r)  =  1(1  -  t),  r  <  2a:  g(r )  =  0,  r  >  2a) 
may  suffice  even  up  to  high  concentrations.  The  the¬ 
oretical  procedure  described  in  this  paper  agrees  very 
well  with  experiments  over  a  wide  range  of  con¬ 
centrations  and  ka  values,  0  <  <■ <,  40"..  and 
0  <  ka  <  10.  Of  course,  these  conclusions  also 
depend  on  the  scattcrcr  properties  relative  to  the  em¬ 
bedding  medium. 
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Acoustic  wave  attenuation  due  to  multiple  scattering  in  a  two-phase  medium  consisting  of  a  fluid 
with  embedded  rigid,  fluid,  or  elastic  particles  of  varying  sizes  is  discussed.  The  formulation, 
involving  the  exciting  and  scattered  fields  of  an  incident  acoustic  plane  wave,  is  based  on  the  T- 
™  matrix  method.  The  propagation  features  of  coherent  waves  in  the  mixture  are  described  by  the 

dispersion  equation  which  is  derived  by  applying  standard  statistical  approximations  to  the 
discrete  random  medium.  Special  attention  is  focused  on  the  pair-correlation  function  between 
the  scatterers  using  the  self-consistent  approximation  (SC A)  which  seems  better  than  the  Percus- 
Yevick  approximation  (PYA)  when  the  volume  fraction  becomes  significant.  Besides  deriving 
low-frequency  analytical  results  for  coherent  wave  speed  and  attenuation,  the  dispersion  equation 
has  been  solved  numerically  for  higher  frequencies  for  particles  with  uniform  and  Gaussian  size 
distributions. 

PACS  numbers:  43.20.Fn,  43.20.  Hq,  43.20.Bi 


INTRODUCTION 

A  study  of  wave  propagation  in  a  multi-component  flu¬ 
id  medium  is  helpful  in  analyzing  the  gross  properties  of  the 
mixture  as  a  whole.  The  propagation  and  attenuation  results 

■  derived  from  a  typical  investigation  yield  data  that  are  valu¬ 
able  in  various  practical  situations. 

In  formulating  the  theory  for  sound  propagation  in  a 
fluid-particle  mixture,  two  fundamentally  different  ap- 
proaches  have  been  followed  in  the  literature.  In  the  first 
approach,  the  solid  particles  are  treated  as  scattering  centers 
and  a  set  of  coupled  equations  is  formulated  to  describe  the 
multiple  scattering.  By  computing  the  acoustic  fields  result¬ 
ing  from  this  scattering  phenomenon,  one  arrives  at  formu¬ 
las  for  the  bulk  parameters  characterizing  the  medium.  In 
the  second  approach,  the  field  variables  characterizing  the 
medium  are  related  through  usual  conservation  laws  and  the 
response  of  the  medium  to  incident  acoustic  pulses  are  de- 
m  rived  from  nonequilibrium  thermodynamic  considerations 
'  •  coupled  with  the  hyperbolic  nature  of  the  governing  system. 
The  results  accrued  from  both  these  approaches  are  largely 
seen  to  complement  each  other  although  a  particular  prob- 
lem  may  entail  selection  of  either  of  the  procedures  in  prefer- 
-*  ence  to  the  other. 

For  suspended  particles,  multiple  scattering  theory 
combined  with  suitable  statistical  approximations  lead  to  re¬ 
sults  that  are  valid  for  a  wide  range  of  frequencies.  This  ap¬ 
proach  takes  account  of  the  microscopic  features  of  the  scat¬ 
terers,  the  effects  of  which  are  finally  reflected  in  the 
coherent  wave  analysis.  Although  rigorous  theoretical  in¬ 
vestigations  of  the  acoustic  wave  attenuation  in  a  fluid-parti¬ 
cle  medium  can  be  traced  to  the  work  of  Sewel,1  systematic 
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multiple  scattering  formulations  capable  of  uniformly  han¬ 
dling  both  the  asymptotic  and  intermediate  frequency 
ranges  and  based  on  realistic  particle  distributions  and  cor¬ 
relation  effects  were  developed  only  in  the  last  few  years. 3-1 1 

In  this  paper  the  multiple  scattering  of  acoustic  waves 
by  suspended  particles  is  described  using  the  T-matrix  (Wa¬ 
terman1'*)  to  characterize  the  single  scatterer  response  and  a 
configurational  average  over  the  random  positions  of  the 
particles.  The  method  presented  leads  to  a  computational 
scheme  that  is  suitable  for  scatterers  of  arbitrary  shape,  ori¬ 
entation,  dense  concentration,  and  at  wavelengths  compara¬ 
ble  to  scatterer  size.  The  complex,  effective  wavenumber  in 
the  random  medium  is  computed  as  a  function  of  frequency. 
It  is  observed  that  the  results  crucially  depend  on  the  volume 
fraction  c,  and  for  c  >  0. 1  the  effects  of  the  pair  correlation 
function  for  the  given  distribution  is  significant  especially  at 
lower  frequencies.  Accordingly,  improved  pair  correlation 
functions  using  the  self-consistent  approximation  (SCA) 
have  been  incorporated  into  the  numerical  algorithm.  The 
SCA  is  a  linear  combination  of  the  Percus-Yevick  (PYA) 
and  the  Hypemetted  Chain  (HNC)  approximation  to  the 
pair  correlation  function.  This  enables  us  to  get  numerical 
results  for  c~0.35.  Closed  form  expressions  for  the  phase 
velocity  and  attenuation  are  presented  in  the  long  wave¬ 
length  limit. 

I.  MULTIPLE  SCATTERING  FORMULATION 

In  the  present  study,  the  scatterers  are  assumed  to  be 
either  rigid,  fluid,  or  elastic  whose  properties  differ  from  the 
embedding  fluid  medium,  which,  for  all  practical  purposes, 
will  be  assumed  here  to  be  inviscid  liquid  akin  to  seawater. 
The  present  formulation  deals  with  the  multiple  scattering 
effects  of  the  solid  phase  on  the  coherent  wave  attenuation. 
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Consider  X  identical,  arbitrarily  shaped  random  scat¬ 
tered  with  smooth  boundary  surface  S  in  an  infinite  fluid 
medium  which  are  referred  to  a  coordinate  system  centered 
at  O  The  points  O,  and  Ol  denote  the  centers  of  the  ith  and 
yth  particles,  respectively;  and  they  are  referred  to  the  origin 
by  the  spherical  polar  coordinates  (r, ,  6, ,  <j>, ).  Pis  any  point  in 
the  medium  outside  the  scattered.  Let  a  time  harmonic 
compressional  wave  of  unit  amplitude  and  frequency  cj  be 
incident  on  the  scattering  medium.  Suppressing  the  time  de¬ 
pendence  of  all  quantities,  we  represent  the  pressure  field 
corresponding  to  the  incident  wave  in  the  form 

6°  =  elkl,  (1) 

where  k  is  the  acoustic  wavenumber,  k  =  oj/cf.  The  sound 
speed  cf  is  given  by  \Ar/pf)' n,  where  A;  and  pf  are  the  com¬ 
pressibility  and  density  of  the  fluid,  respectively.  We  denote 
the  compressional  and  transverse  wave  speeds  in  the  elastic 
scatterer  by  and  c,  given  by  cp  —  [(X  +  2 p)/p\u~  and 
c,  =  [( p/p)]1'-,  respectively,  where  A  and  p  are  Lame's  con¬ 
stants  and  p  is  the  density.  For  a  fluid  scatterer,  cp  =  [(/l  / 
p)]l/2,  where  A  an  dp  are  the  compressibility  and  density  of 
the  fluid,  respectively. 

The  total  pressure  field  at  any  point  outside  the  scat- 
terers  is  given  by, 

<5  |r|  =  <#  °(r)  +  £^;(r-r,|,  (2) 

j  =  i 

where<3y‘(r  —  r;)is  the  field  scattered  by  theyth  particle  to  the 
point  of  observation  r. 

We  now  observe  that  the  field  d>‘  exciting  the  yth  scat¬ 
terer  is  the  resultant  of  the  incident  field  4  0  and  the  scattered 
field  from  all  other  scatterers  so  that 

<*;ir)  =  <*°(r)  +  2<*;(r-r,);  a<|r  -  r,  |<2a,  (3) 

where  a  is  the  radius  of  the  sphere  circumscribing  the  solid 
particle  and  the  superscript  e  refers  to  the  exciting  field.  We 
also  assume  that  the  transparent  spheres  of  radius  a  circum¬ 
scribing  each  particle  do  not  interpenetrate. 

The  multiple  scattering  formulation  developed  here  is 
based  on  the  T-matrix  approach,  see  for  example,  Varadan 
and  Varadan.15  We  first  expand  the  field  quantities  in  terms 
of  spherical  wavefunctions 

[Ou  1  f  h,[kr)\ 

L.H-LK  141 

where  h,  are  the  spherical  Hankel  functions,  Ylma  the  nor¬ 
malized  spherical  harmonics,  and  j,  the  spherical  Bessel 
functions;  a  —  e  or  o  refers  to  the  even  or  odd  parity  of  the 
angular  dependence,  /  =  0,1,...;  m  =  0,1,.../. 

As  usual,  the  field  quantities  that  are  regular  at  the  ori¬ 
gin  will  involve  the  spherical  function  j,  instead  of  h,.  We 
thus  write 

<*.'lrl=  £  2  ][/LOu/'Llr-rJ.  <5| 

i  *  O  m  -  O  tj  ,  • 

d;trl=2  £  2  a‘i„„  Re  (/''m„(r-r/),  16! 

f  i  O  m  -  0  „  -  • 


*>)=£  I  Zg'L,  Ren^(r).  (7| 

I  O  m  0 

where  g°m<7  is  the  known  incident  field  coefficient, /;m„  and 
a,ma  are  the  unknown  scattered  and  exciting  field  coeffi¬ 
cients,  respectively.  The  exciting  and  scattered  fields  are 
then  related  by 

Re  iF'-Jr-r,) 

Ima 

=  *  >1  +  £  I /w  Ou  'w(r-r, >.  (8) 

i  Ima 

The  T  matrix15  can  be  used  to  relate  the  exciting  and 
scattered  field  coefficients  pertaining  to  a  single  scatterer  as 
follows 


f  Ima  & 

/  'm‘& 


The  right-  and  left-hand  sides  of  Eq.  (8)  refer  to  two 
different  origins  (the  centers  of  the  tth  and  yth  scatterers, 
respectively).  This  is  remedied  by  invoking  the  translation- 
addition  theorem  for  the  spherical  wave  functions16  and  the 
orthogonality  of  these  functions  can  be  used  to  extract  a  set 
of  linear  matrix  equations  for  the  exciting  field  coefficients. 
The  details  of  these  steps  have  been  described  in  Refs.  17-19. 
Performing  the  above  operations  and  further  assuming  that 
the  scatterers  are  identical  we  obtain  for  N— +  oo , 

.  ikk^r.c 

^*\ma 


+  (7V_1)£  2  Trm.&,.n.a.a).„.a. 

Ima  /  "m~a“ 

Xo-,.m<,.w(r,  -  r, ).  (10) 

where  erf r,  —  r, )  is  the  translation  matrix,  the  detailed  form 
of  which  is  not  presented  here. 

If  the  scatterers  are  not  identical  but  there  is  a  distribu¬ 
tion  of  sizes,  then  we  may  replace  Tin  Eq.  ( 10)  by  the  average 
T  matrix  ( T  >  defined  as 

(T)  =  ^T(a)q{a)da,  (11) 

where  q(a)  is  a  function  specifying  the  size  distribution  of 
particles  assuming  that  they  are  all  of  the  same  shape.  In 
order  to  further  refine  the  calculations  one  may  want  to  per¬ 
form  the  average  in  Eq.  (11)  only  after  introducing  the  con¬ 
figurational  average  since  the  nearest  distance  between  par¬ 
ticles  is  a  function  of  their  size.  To  a  first  approximation,  we 
can  use  <r>  in  Eq.  (10)  and  use  a  the  mean  radius  to  describe 
the  particle  size. 


II.  CONFIGURATIONAL  AVERAGE  AND  THE  PAIR 
CORRELATION  FUNCTION 

For  a  system  with  a  large  number  of  scatterers,  it  is 
more  meaningful  to  study  the  effective  propagation  charac¬ 
teristics  in  the  medium  rather  than  the  details  of  the  multiple 
scattering  processes  that  take  place.  Thus  a  configurational 
average  is  performed  in  Eq  ||0|  over  the  positions  of  all 
particles  except  theyth  which  is  assumed  to  be  held  fixed.  We 
thus  have 
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1  ‘  <oU>,  -  +  (.V-  1)  X  S  Trm.a;,.m.a. 


I'm  a  l 0 m’a* 


'  j  (a‘fm-om  n 


./-,„( r,  -  r,]p(r,|r,)<fr,, 


where  <  ),  denotes  the  configuration  average  with  the  /th 
scatterer  held  fixed  and  p\r, ,  r, )  is  the  conditional  probability 
distribution  function. 

Equation  (12)  is  a  hierarchy  which  when  iterated  will 
involve  higher  conditional  probability  distribution  func¬ 
tions.  The  hierarchy  is  truncated  by  invoking  the  quasi-crys¬ 
talline  approximation  (QCA)  first  suggested  by  Lax20  which 
works  surprisingly  well  for  a  wide  range  of  concentrations. 
According  to  the  QCA 

)ij  —  (a‘l  ’m’o’  )(•  (13) 

To  study  the  coherent  or  average  field  in  the  effective 
medium,  we  assume  that  the  average  field  is  a  plane  wave 
propagating  in  the  same  direction  k0  —  z  as  the  original 
plane  wave  but  with  a  complex  propagation  constant 
K  —  Kx  +  iKj  which  is  frequency  dependent.  The  real  part 
Af ,  is  related  to  the  phase  velocity  and  AT,,  the  imaginary  part 
is  proportional  to  the  attenuation  constant.  Thus 

<«;„*>,  ~X,mae'Kl\  (14) 

In  order  to  complete  the  integration  in  Eq.  ( 1 2),  the  joint 
probability  function  must  be  specified.  It  is  convenient  to 
write 


g{x)/V,  x>  1 


where  we  have  assumed  that  the  particles  are  impenetrable 
and  that  for  a  translationally  invariant  system,  Pit,  |r()  de¬ 
pends  only  on  j  r,  —  ry  |  =  2 ax  where  la  is  the  hard  core  radi¬ 
us  or  the  minimum  distance  between  particles,  each  of  radius 
a.  In  the  statistical  mechanics  literature  g[x)  is  known  as  the 
radial  distribution  function. 

Equations  (13H15)  are  substituted  in  Eq.  (12)  and  the 
extinction  theorem  can  be  invoked  to  cancel  the  incident 
wave  term  on  the  right-hand  side  of  Eq.  (12)  with  part  of  the 
second  term  (refer  to  Ref.  18  for  details).  The  resulting  equa- 


Xla,  ~  X  t.rofXfot 

X  £  /'  ~'(2/  +  \]a\0,l'\0,l\A  )JHa 


X(!  6c/((Acal2  —  (Afa(2l(  +  24c/A), 


where 


JHX  =  2kaj,\lKa\h\[2ka)  -  2Kahx[2ka)jx\2Ka), 

.  fl(o,/'jo./:/i)  =  i2x  +  n("<]2. 

|_  where  [ioo]  is  the  Wigner  3-y  symbol.16 


(  g(-*l  -  1  \  j ^2Kax)hx\2kax\dx, 


and  c  =  l4tr/}'a\\ /y  is  the  fractional  volume  occupied  by 

the  particles. 


Equation  (16)  is  a  system  of  linear  simultaneous  equa¬ 
tions  for  the  coefficients  X,.  For  a  nontrivial  solution  of  the 
coherent  field  we  must  require  the  determinant  of  the  coeffi¬ 
cient  matrix  to  vanish.  This  is  the  required  dispersion  equa¬ 
tion  which  can  be  solved  for  the  effective  propagation  con¬ 
stant  K  as  a  function  of  =  ca/c/  and  c .  The  determination  of 
K  from  Eq.  (16)  is  necessarily  numerical  except  in  the  long 
wavelength  limit.  This  will  be  discussed  in  the  next  section. 
We  now  proceed  to  discuss  the  forms  of  g(x)  that  will  be  used 
to  numerically  evaluate  Ix . 

Several  models  ofg(x)  are  available.  For  uncorrelated, 
impenetrable  particles 


1/(1  -c), 

0, 


x>  1 
X<  1. 


(20) 


This  approximation  for  g(x)  refers  to  Talbot  and  Willis,21 


known  as  the  well-stirred  approximation  (WSA),  and  is  ex¬ 


pected  to  be  valid  for  low  values  of  c,  and  as  discussed  by 
Bringi  et  al fails  at  c  >  0. 1 25.  For  concentrations  c  <  0. 1 25, 
there  is  an  analytical  form  forg(x)  due  to  Matem.22  Twersky8 


has  used  a  virial  expansion  to  obtain  g(x)  shown  as 


0,  x  <  1 

g(xj  =  •  1  +  8c(l  —  jx  +  3) .  1  <x  <  2  (21) 

\,  x>2 


which  is  valid  at  low  concentrations. 

Improved  models  of  the  pair  correlation  function  valid 
for  concentrations  up  to  40%  are  the  Percus-Yevick  approx¬ 
imation  (PYA)  and  the  Hypemetted-Chain  approximation 
(HNC).  The  Percus-Yevick  model25  has  been  solved  analyti¬ 
cally  by  Wertheim24  for  the  case  of  hard  impenetrable  parti¬ 
cles.  It  is  expected  to  be  somewhat  better  than  the  HNC  25 
One  of  the  defects  of  the  PYA  is  that  the  two  equations  that 
can  be  derived  for  the  pressure  P  in  a  fluid  containing  "hard” 
particles  lead  to  different  answers  when  the  PYA  forg(x)  are 
substituted  in  them  (these  equations  are  derived  by  Percus23). 
Rowlinson:f>  remedied  this  by  assuming  that  the  direct  cor¬ 
relation  function  which  is  the  short  range  part  of  the  correla¬ 
tion  function  is  a  linear  combination  of  the  ones  resulting 
from  the  PYA  and  HNC  models.  They  were  combined  with 
an  adjustable  parameter  d>  and  the  two  pressure  equations 
were  solved  simultaneously  for  P  and  <b.  This  is  called  the 
self-consistent  approximation  (SC  A)  and  it  is  valid  for  higher 
concentrations  than  the  PYA  and  HNC  models. 


III.  RAYLEIGH  LIMIT  SOLUTION 

The  dispersion  relation  derived  from  Eq.  (16)  can  be 
solved  in  detail  to  predict  the  attenuation  features  and  wave 
speeds  of  coherent  acoustic  waves  in  the  two-phase  medium. 
Although  the  system  of  equations  requires  a  numerical  ap¬ 
proach  to  yield  solutions  for  higher  values  of  frequency,  ana¬ 
lytical  results  can  be  obtained  for  low-frequency  approxima¬ 
tions.  Including  the  effects  of  correlation  between  particles  it 
is  seen  that  an  attenuation  factor  is  obtained  even  in  the  low- 
frequency  approximation.  Analytical  results  are  seen  to 
mainly  depend  on  the  form  of  the  correlation  assumed.  Us¬ 
ing,  for  example,  the  g(x)  given  by  the  virial  series,  Eq.  1 2 1 ), 
and  for  leading  order  in  (fail,  dispersion  equations  can  be 
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derived  from  Eq.  (16)  for  all  the  three  types  of  scatterers, 
namely  rigid,  fluid,  and  elastic  spheres.  We  get  for  the  real 
and  imaginary  parts  of  the  effective  wavenumber 

(*,/*  )2  =  ( 1  +  ac)(  1  -  0c)/(  1  +  20c)  (22) 

and 

|  2 (KJKX)  =  (Ara)3(l  -  8c  +  34c2) 

X [a2/(  1  +  ac)  +  20 2/(  1  +  20c){\  -0c)]/ 3, 

(23) 


where 


N 


For  a  uniform  size  distribution  G  (5)  =  1 ,  while  for  the  Gaus¬ 
sian  size  distribution  G  (a)  is  given  by 


0@|-2(,+  MfMl£)/2,,r 


(-sM 

+3(f)e,p(- 

+o-5[,+"<^)]- 


(25) 


where  m  is  the  standard  deviation.  When  the  concentration 
is  higher  than  about  5%  the  term  ( 1  —  8c  +  34c2)  is  replaced 
by  ( 1  —  c)V(  1  +  2c)2  in  Eq.  (23)  to  correct  for  the  higher  con- 
centration  using  the  Omstein-Zernike  equation8  forg(x).  At 
low  frequencies  the  attenuation  is  due  to  multiple  scattering 
alone  since  the  scatterers  are  assumed  lossless.  However, 
/!■  losses  due  to  other  factors  such  as  viscosity,  friction,  etc., 
-*■  may  dominate  the  scattering  losses  (which  are  proportional 
to  k  3a3)  at  low  frequencies.  At  higher  frequencies,  the  multi- 

V  pie  scattering  losses  increase  significantly  and  may  dominate 
-  the  viscous  or  frictional  losses.  Methods  of  introducing  these 

additional  loss  mechanisms  into  the  current  theory  are  un- 
• .  der  investigation  by  the  authors. 

*  IV.  PHASE  VELOCITY  AND  ATTENUATION  AT  HIGHER 
FREQUENCIES 

V  The  analytical  expressions  for  wave  speed  and  attenu¬ 
ation  factor  as  obtained  above  could  be  derived  only  for  very 


Concentration  c 


FIG.  t.  Coherent  wave  attenuation  versus  concentration  for  ka  =  0.05  for 
rigid  spheres  embedded  in  water. 


low  values  of  ka  as  higher  approximations  would  lead  to 
unwieldy  expressions.  A  quantitative  estimate  of  the  multi¬ 
ple  scattering  process  at  resonant  and  higher  frequencies  can 
be  obtained  by  numerically  solving  the  dispersion  equation, 
Eq.  (16).  In  fact,  detailed  numerical  results  can  be  obtained 
for  rigid,  fluid,  or  elastic  scatterers  of  different  geometries  as 
the  only  additional  input  parameter  in  such  cases  is  the  ap¬ 
propriate  T  matrix.  The  computational  scheme  has  been  de¬ 
scribed  previously  and  will  not  be  repeated  here.27  28 

Essentially,  the  solution  of  the  dispersion  equation  de¬ 
rived  from  Eq.  (16)  valid  at  higher  frequencies  involves  an 
iterative  procedure  for  determining  the  dominant  root  in  the 
complex  plane  in  terms  of  ka.  Although  the  algorithm  used 
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FIG.  2.  Attenuation  in  dB/ft  versus  concentration  for  rigid  spheres  for 
ka  =  0.6. 
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FIG.  6.  Relative  phase  velocity  versus  concentration  for  spherical  “red 
blood  cells"  in  a  plasma  medium  for  ka  =  0.6  and  1. 0  using  SCA. 
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FIG.  3.  Attenuation  coefficient  versus  concentration  for  rigid  spheres  in 
water  for  ka  =  I  0. 
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.  4.  Relative  phase  v  elocity  '  ersus  concentration  for  ka  =  0.05  and  1 .0. 
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FIG.  7.  Relative  phase  velocity  versus  concentration  for  rigid  spheres  with 
Gaussian  size  distributions. 
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FIG.  5.  Attenuation  coefficient  versus  concentration  for  spherical  “red 
blood  cells"  in  a  plasma  medium  for  ka  =  06  and  I  0  using  SCA. 
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FIG  8  Attenuation  coefficient  versus  concentration  for  rigid  spheres  with 
Gaussian  size  distributions. 
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appears  to  be  insensitive  to  considerable  deviations  from  the 
true  solution,  good  initial  values  are  provided  by  the  low- 
frequency  results  of  the  previous  section.  As  expected,  the 
order  of  the  coefficient  matrix  that  has  to  be  retained  in  the 
amplitude  equations  for  a  convergent  solution  is  directly 
proportional  to  the  frequency.  The  same  is  also  true  with 
respect  to  the  size  of  T  matrix. 

In  Fig.  1  we  show  the  imaginary  part  of  the  propagation 
constant  K/k  for  rigid  spheres  immersed  in  water  with 
ka  =  0.05  where  k  is  the  wavenumber  of  the  incident  acous¬ 
tic  wave  in  water.  The  curve  depicts  the  correct  functional 
form  for  K2  vs  c  and  shows  a  broad  maximum  in  the  range 
0.15  <c  <0.20.  At  higher  concentrations,  the  attenuation 
decreases  implying  that  the  random  medium  appears  more 
homogeneous  with  "restricted**  randomness.  In  the  limit  full 
packing  (c~ 0.64  for  spheres),  the  attenuation  should  attain  a 
fixed  value  consistent  with  a  bulk  rigid  medium.  The  form 
forgfjc)  in  the  limit  or  full  packing  is  not  known,  hence  com¬ 
putations  could  not  be  performed. 

In  Fig.  2  we  show  a  similar  plot  of  coherent  attenuation 
versus  concentration  for  rigid  spheres  in  water  for  compari¬ 
son  with  measurements  conducted  by  Hampton.24  Note  that 
the  ordinate  scale  is  in  dB  per  foot  where  dB/ft  =  1 109 (AT,/ 
k ).  The  value  of  ka  is  chosen  to  be  0.6  so  as  to  yield  computed 
attenuation  values  in  the  range  observed  by  Hampton  (his 
experimental  points  are  shown  by  circles  in  Fig.  2).  His  ex¬ 
periments  involved  clay  particles  in  a  polydisperse  distribu¬ 
tions  embedded  in  water  with  maximum  ka  values  much 
smaller  than  0.6.  Computation  at  such  small  ka  values  yield 
attenuations  which  are  much  smaller  than  observed  since  AT, 
is  of  O  {k  \jJ).  We  conclude  that  the  experimentally  observed 
attenuation  is  due  to  loss  mechanisms  other  than  multiple 
scattering.  In  any  case  the  present  computations  agree  favor¬ 
ably  with  experimental  values  at  ka  =  0.6  showing  the  ob¬ 
served  functional  dependence  with  increasing  concentra¬ 
tion.  Note  that  the  computed  attenuation  is  due  to  multiple 
scattering  losses  only.  For  c  >  0.05,  the  single  scattering  ap¬ 
proximation  significantly  overestimates  the  attenuation,  in¬ 
creasing  without  limit  as  c  approaches  dense  packing.  The 
self-consistent  form  forg(x)  appears  most  suitable  for  the  full 
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FIG  4  Relative  phase  velocity  versus  concentration  for  elastic  spheres  wiih 
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FIG.  10.  Attenuation  coefficient  versus  concentration  for  elastic  spheres 
with  Gaussian  size  distributions. 


range  of  concentrations.  In  Fig.  3  we  again  show  K2/k  versus 
concentration  for  ka  =  1.0.  Figure  4  shows  the  phase  veloc¬ 
ity  in  the  bulk  medium  relative  to  the  phase  velocity  of  the 
background  medium  in  which  the  scatterers  are  embedded. 
At  higher  concentrations  there  is  evidence  of  dispersion  as 
seen  by  the  separation  of  the  ka  =  0.5  and  ka  =  1.0  curves. 
The  bulk  or  composite  medium  supports  “slow”  acoustic 
waves  when  the  scatterers  are  rigid. 

Computations  were  also  performed  for  fluid  and  elastic 
scatterers  with  the  only  change  in  the  computation  occur¬ 
ring  as  a  result  of  the  appropriate  T  matrix'0  51  in  Eq.  (16). 
The  acoustic  parameters  of  the  fluid  scatterers  and  the  sur¬ 
rounding  fluid  medium  were  chosen  to  bep  =  1 .092  gem-3, 
cp  =  1.64  X  105  cm/s,  and  pf  =0.021,  cf  =  1.55  X  105,  re¬ 
spectively,  where  p  is  the  density  and  c  is  the  sound  veloc¬ 
ity.'2  These  parameters  are  appropriate  for  red  blood  cells 
(RBC)  in  isotonic  plasma  although  a  spherical  model  for  the 
RBC  is  an  approximation.  Figure  5  shows  the  computed 
attenuation  coefficient  as  a  function  of  concentration  for 
ka  =  0.6  and  ka  =  1 .0  using  the  self-consistent  pair-correla¬ 
tion  function.  Note  that  the  computed  attenuation  is  due  to 
multiple  scattering.  In  Fig.  6  we  show  the  relative  phase  ve¬ 
locity  versus  c  for  ka  =  0.6  and  1.0  using  the  self-consistent 
model.  The  bulk  medium  now  appears  nondispersive  and 
supports  ‘‘fast  waves.”  In  Fig.  7,  the  phase  velocity  versus 
concentration  c  is  shown  for  Gaussian  size  distribution  [m/ 
a  =  0. 1 8  and  0.4)  while  the  corresponding  coherent  attenu¬ 
ation  curves  are  shown  in  Fig.  8.  In  Figs.  9  and  10,  the  phase 
velocity  and  coherent  attenuation  for  elastic  scatterers  with 
Gaussian  and  uniform  size  distributions  are  depicted.  The 
properties  of  the  elastic  scatterers  are  taken  as  cp/cK  ~  2.04. 
p/p,  =  2.7. 
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appears  to  be  insensitive  to  considerable  deviations  from  the 
true  solution,  good  initial  values  are  provided  by  the  low- 
frequency  results  of  the  previous  section.  As  expected,  the 
order  of  the  coefficient  matrix  that  has  to  be  retained  in  the 
amplitude  equations  for  a  convergent  solution  is  directly 
proportional  to  the  frequency.  The  same  is  also  true  with 
respect  to  the  size  of  T  matrix. 

In  Fig.  1  we  show  the  imaginary  part  of  the  propagation 
constant  K/k  for  rigid  spheres  immersed  in  water  with 
ka  =  0.0S  where  k  is  the  wavenumber  of  the  incident  acous¬ 
tic  wave  in  water.  The  curve  depicts  the  correct  functional 
form  for  K2  vs  c  and  shows  a  broad  maximum  in  the  range 
0.15<c<0.20.  At  higher  concentrations,  the  attenuation 
decreases  implying  that  the  random  medium  appears  more 
homogeneous  with  “restricted”  randomness.  In  the  limit  full 
packing  (c~0.64  for  spheres),  the  attenuation  should  attain  a 
fixed  value  consistent  with  a  bulk  rigid  medium.  The  form 
forgfoc)  in  the  limit  or  full  packing  is  not  known,  hence  com¬ 
putations  could  not  be  performed. 

In  Fig.  2  we  show  a  similar  plot  of  coherent  attenuation 
versus  concentration  for  rigid  spheres  in  water  for  compari¬ 
son  with  measurements  conducted  by  Hampton.29  Note  that 
the  ordinate  scale  is  in  dB  per  foot  where  dB/ft  =  1 109(AT2/ 
k ).  The  value  of  ka  is  chosen  to  be  0.6  so  as  to  yield  computed 
attenuation  values  in  the  range  observed  by  Hampton  (his 
experimental  points  are  shown  by  circles  in  Fig.  2).  His  ex¬ 
periments  involved  clay  particles  in  a  polydisperse  distribu¬ 
tions  embedded  in  water  with  maximum  ka  values  much 
smaller  than  0.6.  Computation  at  such  small  ka  values  yield 
attenuations  which  are  much  smaller  than  observed  since  K2 
is  of  O  (k  V).  We  conclude  that  the  experimentally  observed 
attenuation  is  due  to  loss  mechanisms  other  than  multiple 
scattering.  In  any  case  the  present  computations  agree  favor¬ 
ably  with  experimental  values  at  ka  —  0.6  showing  the  ob¬ 
served  functional  dependence  with  increasing  concentra¬ 
tion.  Note  that  the  computed  attenuation  is  due  to  multiple 
scattering  losses  only.  For  c  >  0.05,  the  single  scattering  ap¬ 
proximation  significantly  overestimates  the  attenuation,  in¬ 
creasing  without  limit  as  c  approaches  dense  packing.  The 
self-consistent  form  forgfjc)  appears  most  suitable  for  the  full 


FIG  9  Relative  phase  velocity  versus  concentration  for  elastic  spheres  with 
Gaussian  si  it  distributions. 

I  94« 


0  0  0  05  010  015  0  20  0  25  0  30 


Concentration  c 


FIG.  10.  Attenuation  coefficient  versus  concentration  for  elastic  spheres 
with  Gaussian  size  distributions. 


range  of  concentrations.  In  Fig.  3  we  again  show  K2/k  versus 
concentration  for  ka  —  1 .0.  Figure  4  shows  the  phase  veloc¬ 
ity  in  the  bulk  medium  relative  to  the  phase  velocity  of  the 
background  medium  in  which  the  scatterers  are  embedded. 
At  higher  concentrations  there  is  evidence  of  dispersion  as 
seen  by  the  separation  of  the  ka  —  0.5  and  ka  —  1.0  curv  es. 
The  bulk  or  composite  medium  supports  “slow”  acoustic 
waves  when  the  scatterers  are  rigid. 

Computations  were  also  performed  for  fluid  and  elastic 
scatterers  with  the  only  change  in  the  computation  occur¬ 
ring  as  a  result  of  the  appropriate  T  matrix30’1  in  Eq.  (16). 
The  acoustic  parameters  of  the  fluid  scatterers  and  the  sur¬ 
rounding  fluid  medium  were  chosen  to  bep  =  1.092  gem-’, 
cp  =  1.64X  105  cm/s,  and  pr  =  0.021,  cf  =  1.55  x  10s,  re¬ 
spectively,  where  p  is  the  density  and  c  is  the  sound  veloc¬ 
ity.’2  These  parameters  are  appropriate  for  red  blood  cells 
(RBC)  in  isotonic  plasma  although  a  spherical  model  for  the 
RBC  is  an  approximation.  Figure  5  shows  the  computed 
attenuation  coefficient  as  a  function  of  concentration  for 
ka  =  0.6  and  ka  =  1 .0  using  the  self-consistent  pair-correla¬ 
tion  function.  Note  that  the  computed  attenuation  is  due  to 
multiple  scattering.  In  Fig.  6  we  show  the  relative  phase  ve¬ 
locity  versus  c  for  ka  =  0.6  and  1.0  using  the  self-consistent 
model.  The  bulk  medium  now  appears  nondispersive  and 
supports  "fast  waves."  In  Fig.  7,  the  phase  velocity  versus 
concentration  c  is  shown  for  Gaussian  size  distribution  I m/ 
a  =  0. 1 8  and  0.4)  while  the  corresponding  coherent  attenu¬ 
ation  curves  are  shown  in  Fig.  8.  In  Figs.  9  and  10.  the  phase 
velocity  and  coherent  attenuation  for  elastic  scatterers  with 
Gaussian  and  uniform  size  distributions  are  depicted.  The 
properties  of  the  elastic  scatterers  are  taken  as  cp/c,  =  2.04. 
p/p,  =  2.7. 
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1  ’  V.  CONCLUSIONS 

This  paper  has  presented  a  T-matrix  formulation  of 
.  multiple  scattering  of  acoustic  waves  in  a  two-phase  random 
'  medium.  The  dispersion  equation,  which  characterizes  the 
propagation  of  coherent  waves  in  the  effective  medium,  is 
■  derived  by  including  the  effects  of  pair-correlation.  Explicit 
•  expressions  for  the  attenuation  coefficient  in  the  Rayleigh 
limit  are  derived  using  the  virial  series  expansion  of  the  radi- 
■’ .  al  distribution  function.  Computations  are  performed  for 
\  rigid,  fluid,  and  elastic  spheres  embedded  in  a  fluid  back¬ 
ground  medium  using  a  number  of  models  for  the  radial 
gl  distribution  function.  At  high  concentrations  the  self-consis¬ 
tent  model  appears  to  give  superior  results  based  on  a  com- 
'•  parison  with  experimental  observations.  The  theory  and 
computational  procedure  developed  in  this  paper  should 
prove  useful  in  dealing  with  acoustic  wave  propagation 
through  ocean-bottom  sediments,  mixtures,  and  suspen¬ 
sions  and  composite  biological  tissues. 
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ABSTRACT 

Acoustic  waves  propagating  through  suspended  materials  with  a  con¬ 
siderable  concentration  (1%<w<40%)  in  the  deep  ocean  environment  is 
examined.  The  particles  in  suspension  have  a  size  distribution  and  their 
relative  positions  are  described  by  the  pair-correlation  function. 
Different  equations  representing  the  pair-correlation  function  are 
investigated.  The  characteristics  of  the  multiple  scattered  waves  are 
presented  as  the  dispersion  of  the  phase  velocity  specified  by  the  real 
part  of  the  effective  wave  number  K  and  the  attenuation  in  wave  intensity 
shown  as  the  loss  tangent. 

The  dependence  of  the  wave  propagation  on  the  size  distribution  and 
concentration  of  particles  is  also  discussed  in  the  low  frequency  range. 
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INTRODUCTION 

An  acoustic  technique  has  been  proposed  to  investigate  the  wave 

propagation  through  the  suspended  sedimentary  particles  in  a  water 

mass^.  The  attenuation  of  acoustic  waves  by  suspended  materials  in 

2-4 

a  fluid  medium  has  been  investigated  extensively  and  the  experi¬ 
mental  attenuation  coefficients  were  found  to  be  compatible  with  those 
derived  from  theory  for  low  concentration  cases.  However,  when  the 
concentration  becomes  higher  the  acoustic  transmission  and  reflection 
will  be  affected  differently  due  to  multiple  scattering.  The  single 
scattering  theory^  which  is  suitable  for  a  sparse  distribution  of  scat- 
terers  (suspended  particles)  is  no  longer  valid  for  a  dense  distribution 
of  scatterers  since  the  higher  order  statistics  known  as  the  pair- 
correlation6  will  be  required  in  the  multiple  scattering  analysis.  In 
other  words,  in  a  medium  containing  a  large  number  of  scatterers  the 
position  of  one  scatterer  is  constrained  by  the  other  therefore  the 
scattered  field  will  be  affected  by  such  a  crowding  characterized  by  a 
correlation  function  among  scatterers. 

The  present  study  examines  acoustic  wave  propagation  through  a 
considerable  concentration  of  sedimentary  particles  suspended  in  water. 

The  quasicrystalline7  assumption  will  be  used  to  truncate  the  hierarchy 
equations  (Foldy-Lax  hierarchy^)  so  that  only  the  pair-correlation  between 
two  particles  is  considered.  As  can  be  seen  the  pair-correlation  function 

g 

is  essential  in  solving  the  coherent  field.  Therefore,  the  hole  correction  , 
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the  virial  expansion  and  the  Percus-Yevick  equations  representing  the 
pair-correlation  function  are  discussed.  However,  it  can  be  shown  that 


all  those  equations  will  yield  the  similar  results  for  moderate  concen¬ 
trations  of  scatterers  in  the  low  frequency  range. 

For  practical  purposes,  a  simpler  form  of  the  pair-correlation 
function  for  elastic  scatterers  is  given  instead  of  the  complex 
Percus-Yevick  equation.  The  scatterers  with  and  without  size  distri¬ 
butions  are  both  treated  in  this  study  to  judge  their  effects  on  the 
wave  propagation  characteristics.  The  results  presented  as  the  dispersion 
of  the  phase  velocity  as  well  as  the  attenuation  against  the  concentration 
are  particularly  in  the  low  frequency  range. 


STATEMENT  OF  THE  PROBLEM 


Consider  a  slab  containing  a  large  number  of  scatterers.  A  time- 
harmonic  plane  wave  of  sound  is  normally  incident  upon  the  slab  (see 
Figure  1)  along  the  negative  z  direction  with  the  time  dependence  e+llvt. 
The  incident  wave  with  a  unit  amplitude  can  thus  be  expressed  as 

t(ft)  -  eik2 

where  k  =  wave  number. 

The  scatterers  are  suspended  particles  modelled  as  elastic  spheres 
whose  properties  ire  shown  in  Table  1.  The  surrounding  acoustic  medium 
(covering  slab)  is  water  with  the  properties  also  shown  in  Table  1.  In 
the  absence  of  scatterers  from  the  medium,  will  then  satisfy  the 

wave  equation 

2  2 

V  'p  +  k  =  0 

However,  the  average  field  (coherent  field)  cijo11  for  any  wave  in  a 
medium  containing  random  particles  can  generally  be  expressed  by 
assuming  it  satisfies  the  following  wave  equation  specified  by  an  ef¬ 
fective  wave  number^  K (sometimes  called  propagation  constant  or  bulk 
parameter),  i.e., 

2  2 

7  <’l>>  +  K  <y>  =  0 

For  low  concentrations  (the  number  of  scatterers  in  the  medium  is 
small  and  the  distances  between  scatterers  are  large  compared  to  the  in¬ 
cident  wave  length)  the  effective  wave  number  K  is  well  known  as1*",  hen 
including  the  size  distribution. 
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K2  =  [k+^-  (f(a,ir)  +  f(a,o))]  [k+^-  (ffa.m)  -  f(a,o))],  (2) 

where  p  =  number  density  (number  of  scatterers  per  volume) 
f(a,ir)  =  J"  f(a,ir)  q(a)da 


f(a,o)  =  fj  f (a,o)  q(a)da 


q(a)  =  the  size  distribution  function 
f(a,ir)  =  the  forward  scattering  function* 
f(a,o)  =  the  backscattering  function* 

Lloyd  and  Berry*3  modified  the  above  equation  using  the  multiple  scattering 
treatment  as  follows 

2  2 

K2  =  k2+47rpf  (a,7r)  +  [-f2(a,Tr)+f2(a,o)  -f*  — ^ f 2  (0 )  de  ] , 

k  0  sin-i  3  aH 

(3) 


where 


f(9) 


CO 


l 


(2n+l) i(-l)n 
1+iC 

n 


Pn(cos  0) . 


.Another  frequently  used  equation  for  the  effective  wave  number  K 

<5 

is  given  as, 

K2  =  k2+4iTof  (a,7r)  +  (47rpf  (a, n) ) 2  fQ  eikR  g(R)dR,  (4) 

where  g(R)  is  the  radial  distribution  function.  Results  using  Eqs. 

(2),  (3)  and  (4)  will  be  examined  in  this  paper. 
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Equation  (2)  can  be  further  reduced  to  the  following  form  provided 


p«l. 


+  4npf(a,ir)  +  0  (p^) 


in  which  0  (  )  represents  the  order  of  (  ).  For  small  scatterers 

and  low  concentrations  a  common  form  used  by  many  investigators  is 
obtained  by  expanding  the  square  root  of  Eq.  (5)  as 


K  =  k  +  2-rrpf  (a,ir)/k 


One  sees  from  Eq.  (6)  that  K  is  actually  a  complex  number  due  to 
the  fact  that  f (a, it)  is  complex.  The  real  and  imaginary  parts  of  K  are 
found  to  be 


Kj^  =  k  +  2TrpfR(a,Tr)/k 


where  the  subscripts  R  and  I  denote  real  and  imaginary  respectively. 
By  using  the  forward  scattering  theorem  which  is 


the  imagninary  part  of  K  (also  known  as  the  attenuation  constant)  can 
thus  be  related  to  the  total  scattering  cross  section  as  as 


2K  =  pa  +  (pa  ) 
is  a 


In  Eq.  (9)  a  is  the  absorption  cross  section  only  if  the  scatterers  are 

2L 

absorbent  (the  absorption  mechanisms  can  be  introduced  by  using  complex 
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internal  wave  numbers  in  the  scattering  functions  f(a,Tr)  and  f(a,o)  or 
it  will  not  appear.  The  absorption  effect  for  sedimentary  particles  is 
small  and  can  be  neglected.  ^ 

Strictly  speaking,  Eqs.  (5),  (6),  (7)  and  (9)  can  only  be  applied 
to  the  low  concentration  cases.  Each  of  the  sparsely  distributed  scatterers 
under  such  circumstance  can  be  treated  as  independent  scatterers,  i.e., 
uncorrelated  and  without  interference.  Therefore  the  total  scattering 
is  just  the  sum  of  the  scattering  from  each  scatterer  which  is  the  basis 
of  the  single  scattering  theory.  A  question  will  arise  as  to  how  large  a 
concentration  will  make  the  single  scattering  invalid?  The  answer 
depends  on  the  type  of  scatterers  being  treated.'*'  However,  a  recent 
experimental  study  by  Ishimaru  et.al.*'0  shows  that  when  the  concentration 
is  greater  than  about  0.1%,  the  experimental  attenuation  constant 
departs  markedly  from  that  calculated  by  single  scattering  theory. 

The  term  concentration  used  in  this  paper  is  defined  as  the 
volumetric  percentage  which  is  the  ratio  between  the  total  volume  of 
scatterers  and  the  volume  of  the  embedding  medium.  Since  the  spherical 
scatterers  are  modelled  as  the  suspended  particles  the  concentration  w  is 
thus  represented  as 

w  =  4-pa/3, 

where  a  is  the  radius  of  the  sphere  and  the  overbar  denotes  the  average 
quantity. 

To  examine  the  effective  wave  number  K  for  a  dense  distribution  of 
scatterers  (high  concentration)  one  needs  to  recall  the  Foldy-Lax  hierarchy 
for  the  coherent  field  holding  more  than  one  scatterers  fixed.  Fortunately, 
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the  problem  can  be  simplified  without  going  into  higher  order  statistics 
by  using  the  quasicrystalline  approximation  to  truncate  the  equations  in 
order  that  only  the  pair  statistics  between  two  scatterers  is  required  in 
evaluating  the  coherent  field.  Such  a  treatment  is  useful  in  which  an 
explicit  approximation  for  the  pair  statistics,  expressed  as  the  pair 
distribution  function  g(R)  can  be  introduced  to  integrate  over  the  total 
volume  to  obtain  the  exciting  coherent  field.  The  second  equation 
in  the  Foldy-Lax  hierarchy  using  the  quasicrystalline  approximation  can 
generally  be  written  as 


=  *.nc(l)  +  p  f  g$rl)f(a,e)  <i|;1(^i)>i  d$.  E(l-i)  (10) 

where  <i^1(^)>_.  =  the  average  exciting  field  of  the  jth  scatterers  with 
its  position  held  fixed 
’</ .  (I$.)  =  the  incident  field  at  R. 

yincv  jJ  j 

t'(a,9)  =  f (a, 9)  q (a)da 

)>^  =  the  average  exciting  field  of  the  ith  scatterer  with 

its  position  held  fixed 

=  propagation  function  of  the  scattered  waves 
d^  =  4-tR^dR  for  the  radially  symmetric  scatterers 


SOLOTION  FOR  RADIALLY  SYMMETRIC  SCATTERERS 


The  effective  wave  number  K  for  a  dense  distribution  of  pair-correlated 
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scatterers  can  be  obtained  using  Eq.  (10)  as  follows 

2n 


K2  =  k2  ♦ 


where 


l  An' 
u  n 


n=o 


n  =  K/k 

CO 

A  =  kf  (a,*)  [l+T  iA  (n)"n+mH  ] 
n  nv  *  L  L  n  nmJ 


m=o 


fn Ca,7t)  =  fn(a,ir)  q(a)da  =  fn 
( 2n+l)i 


f.n(a-')  =  1+iC 


n 


n+m| 


.  .  m+n  m 

H  .'l  ■  d.  (°;°)  [i™i  *  X  ) 

nm  L  n  m  L  ,  3  2  . 

'  ■  k  n  -1 

:  n-m ; 


nr 


N  =  4ttP  f  *  [g(R)-l]j  (KR)h2  (kR)R2dR 
m  mm 


m  m 


V1  .  v2 
=  N  +  l  N 
m  m 


j m ( KR )  =  spherical  Bessel  function 


(ID 


(12) 


h“(kR)  =  spherical  Hankel  function  of  the  second  kind 
m 


g (R)  =  pair-correlation  function 


,o.o. 


dm(n;m)  =  Lerjendre  expansion  for 


n+m 


P  (x)P  (x)  =  V  dq(°;°)P  (x) 
n  m  ^  n  m  q 


q=o 


P  (x),  P  (x),  P  (x)  =  Lengendre  polynomials 
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For  the  low  frequency  case,  only  the  monopole  (each  scatterer,  i.e. 
sphere,  is  compressed  and  expanded  by  the  incident  condensations  and  rare 
factions  and  a  spherical  wave  is  thus  radiated  -  this  monopole  type 
radiation  is  dominated  by  the  compressibility  of  the  scatterer  and  is 
independent  of  directions)  and  dipole  (the  scatterer' s  inertia  causes  the 
scatterers  to  have  a  motion  which  is  equivalent  to  the  surrounding  medium 
being  at  rest  and  the  scatterer  being  in  oscillation  -  this  dipole  type 
reradiation  is  affected  by  the  density  of  the  scatterer  and  depends  on 
the  scattering  angle  8)  terms  are  important  in  doing  the  computation. 
Therefore  only  two  terms  in  Eq.  (11)  are  required  and  a  closed  form  solution 
can  be  obtained.  However,  more  terms  are  necessary  in  the  high  frequency 
range  which  make  the  calculation  for  K  a  tedious  job  and  only  a  numerical 
approach  is  possible  at  the  present  time.  This  can  be  done  by  the 
computer  for  a  wide  frequency  range  by  selecting  a  suitable  number  terms 
to  achieve  convergence.  The  acoustic  scattering  analysis  for  suspended 
particles  in  the  high  frequency  range  will  be  discussed  in  a  later  paper. 
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Now  one  can  compare  Eq.  (11)  with  Eq.  (5)  to  see  the  difference.  To 
study  the  wave  propagation  through  a  sparse  distribution  of  scatterers 
does  not  involve  the  iterative  solution  for  K  which  is  instead  required 
for  a  dense  distribution  of  pair-correlated  scatterers.  The  effective 
wave  number  K  in  a  medium  containing  scatterers  of  low  concentrations  is 
essentially  a  slight  perturbation  of  Eq.  (11)  with  n=l.  One  also  sees 
from  Eq.  (12)  that  in  order  to  obtain  Nm  the  pair-correlation  function 
(or  the  radial  distribution  function  for  the  radially  symmetric  scatterers) 
g(R)  needs  to  be  specified  first.  If  there  is  a  closed  form  expression 
for  g(R),  the  computation  can  be  made  much  easier  depending  on  how  simple 
the  form  is. 


11 


RADIAL  DISTRIBUTION  FUNCTION 

Several  radial  distribution  functions  (for  the  radial  symmetric 

scatterers  the  pair  correlation  function  g(R)  is  called  the  radial 

14 

distribution  function  g(R))  have  been  applied  successfully  to  analyze 

the  multiple  scattering  of  waves  by  a  random  distribution  of  pair-correlated 

scatterers.  The  commonly  used  hole  correction  equation  has  been  proved 

9 

to  be  a  poor  approximation  for  appreciable  concentrations.  Twersky  used 
the  virial  expansion  equation  which  is  an  iterative  solution  for  the 
Percus-Yevick  equation  for  impenetrable  scatterers  to  obtain  the  effective 
wave  number  K  for  the  low  frequency  range  and  moderate  concentrations. 

To  make  the  calculation  easier  the  following  form  is  proposed  to 
avoid  a  more  complex  form  in  the  expression  of  the  Percus-Yevick  equation. 
That  is 

eQ  >  R  <b' 

"ll+A'Me  B’(w)(R"b,)cosC’(w)(R-b')  ,  R  >  b'  (13) 

where  A'(w),  B'(w)  and  C'(w)  are  functions  of  the  concentration  w  and  b' 
the  separation  distance  (exclusion  length)  between  two  scatterers.  For 
different  scatterers  with  different  distributions  A'(w),  B'(w)  and  C'(w) 
have  different  forms  and  are  decided  by  the  distribution  measurements 
from  the  field  data.  Generally  speaking,  the  magnitudes  of  A'(w)  and 
C' (w)  increases  with  the  increasing  concentration  and  that  of  B'(w) 
instead  of  increasing  decreases  as  the  concentration  increases. 

The  idea  in  establishing  Eq.  (13)  for  g(R)  is  from  the  experimental 
histogram  of  the  radial  distribution  function  for  the  manganese  nodule 
fields1^  (nodules  are  elastic  materials)  and  the  Percus-Yevick  equation. 


■■r.  vj  .-*  v.  r; 


.y 
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One  sees  that  g(R)  has  the  trend  of  decaying  and  oscillating  about  the 
value  one  when  the  radial  distance  R  is  greater  than  the  exclusion  length 
b(=2a).  The  asymptotic  value  of  g(R)  is  one  when  R  is  generally  three 
or  four  times  larger  than  the  exclusion  length  (see  Figure  2). 

For  a  moderate  concentration,  e.g.  w»0.1,  a  comparison  of  g(R)  using 
different  expressions  is  presented  as  Figure  3.  As  can  be  seen  from  the 
figure,  (in  the  range  of  0<R<2b)  the  curves  representing  the  Percus-Yevick 
equation,  the  virial  expansion  and  the  proposed  radial  distribution  function 
which  has  the  following  form  are  very  close  to  each  other. 


f  0  ,  R< 

^  l+0.3e-7r(R"b 


b 

b)/b 


(for  w  =  0.1) 


g(R)  = 

I  -  CD  K1  /l, 

cos  Tr(R-b)/2b  ,  R  >  b 
As  a  matter  of  fact,  the  virial  expansion  of  g(R)  is  the  simplest  one  in 
doing  the  calculation  with  a  moderate  concentration. 

To  calculate  the  effective  wave  number  K,  the  equation  of  Nm(Eq.(12)) 
which  involves  the  radial  distribution  function  g(R)  needs  to  be  obtained 
first.  In  general,  Nm  can  be  expressed  as  the  following  integral. 


N  =  47rpF(K,k)  f  °°  [g(R)-l]RadR 

m  J 


(14) 


where  F(K,k)  is  the  function  of  the  effective  wave  number  K  and  the  incident 
wave  number  k.  Eq.  (14)  is  generated  due  to  the  ascending  series  expansion 
of  jm (KR)  and  h^  (kR)  in  Eq.  (12)  as  follows 

.  m  _  yn  Jn  r  (-l)v(m+v)  l  X2v 

JmW  “  1  v!  (2m+2v+l)  f 

v=o 


hm  W  -  VX)  -  1  ym(X) 


2 


(15) 
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2v 


y  _  _ L_  V  r (2m-2v+l)X 

nr  J  ,,mvm+l  L  m!r(ra-\>+l) 

2  X  v_ o 

where  r  (  )  is  the  Gamma  function. 

In  general,  the  integral 

I®  ■/"  [g(R)-l]RadR 
J  o 

can  be  expressed  as  a  recurssion  formular  which  is  useful  in  obtaining 
the  iterative  solution  using  the  computer.  A  recurssion  formular  of 
Eq.  (16)  using  Eq.  (13)  for  g(R)  can  be  expressed  as 

/.  i  >  1 

1°  =  ia-  iiil 

a+1  , 

where 


Ct  f a 

T  =/ 

•'V.' 


“  e'B'CR"b,)cos  C'(R-b')RadR 


— — j  [B'(b')01  +  B’aTa"1-C,aT*(ct"1)] 
B'  +C’ 


T*(a-1)  e-B*(R-b')sinC,(R.b,)Ra-l 

*  h  * 


dR 


T°  =  -B' 


.2*r.2 


B'  +C 


-C' 


B,Z+C'2 

For  moderate  concentrations  the  expression  of  Ia  using  the  virial 
expansion  for  g(R)  is 

.a  ,  b3*1  .  r(2b)a+1  3b‘,*1(I-2°‘*2)  ba+1  (2a+4-l) , 


(16) 


where  b  =  2a.  For  the  lowest  order  of  m  in  Eq.  (15),  which  is  useful  in 
the  low  frequency  scattering  for  monopole  and  dipole  case,  N*  is  obtained 
as 

n1  =  4*pi2  (; 

o 

2 

=  - 8w  +  34w 

which  is  the  same  result  as  Eq.  (28)  of  Ref.  9. 

For  the  low  frequency  scattering,  an  exact  form  of  by  using  the 
Ornstein-Zerwike  equation  is 

n‘  =  4,pj T”  [g(R,-l]R2dR  =  -1  C 

which  can  be  shown  to  be  the  same  as  Eq.  (17)  when  the  higher  order  terms 
(0(w3))  are  neglected  after  the  expansion  of  Eq.  (18). 
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APPROXIMATION  IN  THE  LOW  FREQUENCY  RANGE 


In  the  present  study  the  suspended  sedimentary  particles  are  in 
general  2  to  4  microns  in  radius  which  make  the  nondimen sional  frequency 
ka  in  the  range  of  0.01  to  0.1  when  using  the  ultrasonic  wave  of  the 
frequency  1.2  to  12  megahertz.  It  is  known  that  in  the  low  frequency 
Rayleigh  limit  (ka<<l)  only  the  monopole  (n=o)  and  dipole  (n=l)  terms 
dominate  the  solution  for  the  effective  wave  number  K  in  Eq.  (11)  .  Thus 
Eq.  (11)  can  be  rewritten  as 

K  =  k  +  — (Aq  +  Ajfi  )  ( 


1  +  a 
k3  0 
1  .  ilfi.  A 
kJ  1 


where 


Ao  ■  kfo  ri*iAoHoo*1W 


Aj  *  kfj  [1  *  iV  H10  *  iAlHll] 


Hoo  ■  NJ  *  lNo2 


hoi  *  Hio  *  Ni  *  iNi2 


H11  *  T  I1?!  i+  V  *  !  'N2  *  iN22] 

k 
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In  the  low  frequency  range,  AQ  and  A^  can  be  approximated  by" 


l-ia0W 


A!  -  ( - — : — -3 

!  +  IIP  al  -  ai.  W 

k3  y  r1 


where 


1  *  4.p J-  [gfRJ-lJl 


ao  '/o’  TTIc^-  <l(a)d! 


a! '/"  mtr  I(a)da 

o  1 

The  limiting  values  of  aQ  and  Ej^  for  ka<<l  are,  respectively 


ao  '  Da0fc(a,i 


fluid  spheres 


a^  =  -SDa^GCEji 


Eq  =  D(-l)G(a)i 


rigid  spheres 


a^  =  -3D  — ^  G(E)i 


a0  =  Da0eG(a)  i 


E^  =  -3Da^G(a)  i 


elastic  spheres 
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where 


D  =  (ka)  /3 

G(a)  =  the  size  distribution  factor 

a0f  =  1/e  -  1 

alf  =  /  (1+2§) 


1 


Oe 


4 

e‘3  el 


-  1 


e  «  g(CL/C0)' 


el 


a  =  mean  size  (radius) 

g  =  density  ratio  between  the  scatterer  and  the  surrounding  medium 
Cq,  Cl,  =  wave  speed  in  the  surrounding  medium,  compressional  wave 
speed  in  the  scatterer,  shear  wave  speed  in  the  scatterer. 

One  sees  from  Eqs.  (19)  and  (20)  that  the  effective  wave  number  K  is  a 
complex  number  and  can  be  written  as 


K  = 


=  1^  +  iKj,  (Kj  >  0) 


In  general,  the  real  part  of  K(K  )  is  much  greater  than  the  imaginary  part 

K 

of  K(Kj)  and  Eq.  (19)  can  thus  be  rewritten  as,  after  separating  the  real 


and  imaginary  part  of  K, 

A..2  .  (1+ws)  (1-wt) 

J  "  l+2wt 


(21) 
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2K  2  2 

-j£  =  Dw(l-Sv,W)  t^-  <■  ^2w^cl,w0] 


where 


fluid 

s=  C-l  )  G(a)  ,  rigid  sphere 


elastic 


fluid 


t  =  (  j)  G  (a)  ,  rigid  sphere 

a. elastic 


For  a  uniform  size  distribution,  i.e.,  G(a)=l,  the  results  of  (21) 

9 

and  (22)  are  identical  to  Eq.  (73)  obtained  by  Twersky. 

Based  on  Eqs.  (21)  and  (22)  the  normalized  phase  velocity  KR/k  and 

the  loss  tangent  2Kj/Kr  are  calculated  for  waves  propagating  through  the 

water  containing  suspended  particles,  modelled  as  elastic  spheres. 

2 

When  the  concentration  is  higher  than  about  5%  the  term  (l-Sw+34w  )  is 
4  2 

replaced  by  (1-w)  /(l+2w)  in  Eq.  (22)  to  correct  for  the  higher  concen¬ 
tration  in  calculating  the  loss  tangent.  The  rigid  scatterers  are  also 
used  in  the  calculation  for  comparison  purpose  and  in  this  case  Eq.  (11) 
becomes 


c2  .  k2  ♦  *  k2  .  ifei  •  I  •  (7.2w)k: 


( l+2w) ' 


which  is  Twersky's  Eq.  (75)  for  uniform  size  scatterers.  The  size 
distribution  factor  is  considered  to  be  from  the  Gaussian  size  distribution 
and  in  this  case  is* 

c<s>  -  ~A  4)2  '">3  *  It2-^1)  a2  *  4=-*  4c|)2(-)  C4) 

'2  *r  a  av2-r  a 


where  m  =  the  standard  deviation. 


+  y(l  +  erf  (-^— )) 
“  v  2m 
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RESULTS  AND  DISCUSSIONS 

The  acoustic  wave  propagating  through  a  fluid  medium  containing 
a  random  distribution  of  suspended  particulats  with  a  considerable 
concentration  is  usually  dispersed  and  attenuated.  The  effective  wave 
number  K  for  densely  distributed  scatterers  can  be  shown  to  be  different 
from  that  for  sparsely  distributed  scatterers  due  to  the  pair-correlation 
between  scatterers.  The  dispersion  of  the  phase  velocity  is  characterized 
by  the  real  part  of  K(K^)  and  the  attenuation  of  the  wave  intensity  can 
be  described  by  the  imaginary  part  of  K(Kj).  The  effective  wave  number  K 
is  frequency  dependent  except  at  the  low  frequency  range  in  which  the 
normalized  phase  velocity  K  /k  is  independent  of  frequencies  as  shown  in 
Eq.  (21). 

The  normalized  phase  velocity  K^/k  decreases  as  the  concentration 
increases  at  the  low  frequency  range  which  can  be  seen  from  Figure  4. 

The  dispersion  of  the  phase  velocity  has  a  dependence  on  the  size  dis¬ 
tribution  also.  One  sees  from  the  plot  that  particles  with  a  higher 
standard  deviation  in  the  Gaussian  size  distribution  make  the  wave 
disperse  faster.  At  the  zero  concentration  there  should  be  no  dispersion. 
Therefore  the  ratio  between  the  effective  wave  number  and  the  incident 
wave  number  becomes  one  as  verified  by  Eq.  (21). 

The  dispersion  of  the  wave  velocity  in  a  medium  containing  the  uniform 
size  scatterers  is  generally  not  so  strong  as  that  in  a  medium  containing 
size  distributed  scatterers.  One  interesting  phenomenon  can  be  seen  from 
Figure  4  that  for  the  uniform  size  rigid  spheres  the  normalized  phase 
velocity  instead  of  decreasi  ,  increases  with  increasing  concentration. 

This  fact  may  be  explained  by  the  high  impedance  of  the  rigid  body. 
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Figure  5  shows  the  relationship  between  the  loss  tangent  2Kj/Kr  and 

the  concentrations  at  ka  =  0.1,  i.e.  low  frequency  range.  As  can  be  seen 

from  the  figure,  the  attenuation  increases  with  the  increasing  concentration 

until  it  reaches  the  maximum  at  the  concentration  of  fifteen  percent  then 

decreases  with  the  increasing  concentration  up  to  forty  percent.  As  a 

matter  of  fact,  this  phenomenon  and  the  saturation  concentration  of  15%  is 

valid  for  all  cases  at  the  low  frequency  range  and  has  also  been  veri- 

13 

fied  by  experiments.  For  scatterers  with  the  Gaussian  size  distributions, 
one  sees  that  the  more  the  particle  size  deviates  from  its  mean  radius, 
the  larger  the  attenuation  is.  Therefore  the  attenuation  of  suspended 
particles  with  the  uniform  size  distribution,  which  is  the  special  case  of 
zero  standard  deviation  in  Eq.  (26),  is  always  smaller  than  that  of 
particles  with  the  Gaussian  size  distribution.  It  can  also  be  noted  that 
the  attenuation  of  fixed  rigid  spheres  is  larger  than  the  present  sedimen¬ 
tary  particles  which  are  modelled  as  elastic  spheres.  One  sees  from  Eq.  (20) 
that  the  shear  wave  of  the  elastic  material  may  affect  the  attenuation  to 
a  certain  degree.  However,  this  effect  does  not  appear  in  the  present 
problem  due  to  its  small  magnitude. 

The  loss  tangent  is  frequency  dependent  and  its  dependence  in  the  low 
frequency  range  for  a  fixed  concentration  (10%)  is  shown  in  Figure  6.  On 
the  log-log  plot  one  sees  that  all  straight  lines  have  the  slope  three  as 
predicted  by  the  Rayleigh  scattering  theory.  The  attenuation  increases 
with  the  increasing  frequency  in  the  Rayleigh  region  (ka«l).  The  scatterers 
with  a  larger  standard  deviation  in  the  Gaussian  size  distribution 
will  reach  the  assigned  attenuation  at  a  lower  frequency.  The  dispersion 
and  the  attenuation  of  waves  using  Equation  (7)  based  on  the  single 
scattering  theory  and  Eq.  (3)  from  Lloyd  and  Berry^  are  also 
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shown  in  Figures  4  and  7  respectively.  One  sees  that  unless  the  concentration 
is  very  low  the  single  scattering  theory  fails  to  predict  the  propagation 
characteristics. 

Lloyd's  equation  predicts  the  wave  dispersion  well  but  fails  in 
determining  the  attenuation.  Eq.  (4)  using  the  hole 

correction  is  badly  affected  by  the  acoustic  properties  of  scatterers. 

This  can  be  seen  from  the  attenuation  for  rigid  and  elastic  scatterers 
(predicts  better  for  rigid  spheres).  However,  it  is  still  not  possible  to 
obtain  the  right  attenuation  unless  the  concentration  is  very  low.  If 
attention  is  paid,  one  can  see  that  there  is  no  correction  for  the  acoustic 
properties  due  to  the  concentration  in  Figs.  (2)  and  (4)  as  it  appears 
implicitly  in  Eq.  (19) (finally  appeared  in  Eqs.  (21)  and  (22))  due  to  the 
consideration  of  the  pair-correlation  between  scatterers.  Besides  the  lack 
of  the  correction  for  a  higher  concentration.  This  may  be  the  reason  why 
the  single  scattering  theory  and  its  modified  equations  predict  a  lower 
dispersion  and  a  much  higher  attenuation. 

The  hole  correction  equation  representing  the  pair-correlation  function 
is  not  suitable  for  the  concentration  greater  than  about  five  percent  as 
can  be  seen  from  Figure  7.  Besides  this  the  hole  correction  equation  will 
generate  negative  attenuation  coefficients  which  are  nonphysical  results 
even  starting  at  a  comparatively  low  concentration. 

When  one  is  interested  in  the  wave  propagation  through  a  higher 
concentration  of  scatterers  in  the  high  frequency  range,  an  iterative 
solution  can  be  obtained  using  a  computer  using  suitable  numerical  methods. 

The  proposed  radial  distribution  function  (Eq.  13))  is  useful  for  pair- 
correlated  elastic  spheres  whenever  the  concentration  dependent  parameters 
are  decided.  However,  this  will  be  left  to  the  next  paper. 
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TABLE  I.  Elastic  Properties  of  the  Sediments  of  the 
Abyssal  Hill  used  in  Calculation. 

(Data  from  Stoll's  paper) 


Range  of  Grain  Sizes 

(a) 

2-4 

pm 

Grain  Density  (p£) 

2.65 

g/cm3 

Shear  Wave  Velocity 

(cT) 

210 

m/sec 

Compressional  Wave  V 

elocity  (C  ) 

L 

3690 

m/sec 

water  density  1  g/cm3 

and  sound  speed 

in  water  1500  m/sec 

are  used. 

m 
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Figure  2.  The  radial  distribution  function  g(R)  from  Percus -Yevick 
Equation  and  Nodule  Histogram 
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Figure  3.  A  comparison  of  g(RJ  using  different  expressions 
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Figure  4.  Phase  velocity  vs.  w  for  Gaussian  and  Uniform  Size 
Distributions 
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Figure  5.  Coherent  Attenuation  vs.  w  for  Gaussian  and  Uniform 
Size  Distributions 


Figure  7,  Coherent  Attenuation  vs.  w  --  Comparison  of  Various  Theories 


APPLICATION  OF  TWERSKY'S  MULTIPLE  SCATTERING  FORMALISM 
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by 

Y.  Mat*  V.K.  Varadarf  and  V.V.  Varadan* 

Wave  Propagation  Group 
Department  of  Engineering  Mechanics 
The  Ohio  State  University,  Columbus,  Ohio  43210 


ABSTRACT 

Acoustic  wave  propagation  through  a  dense  suspension  of  solid  elastic 
particles  in  water  is  studied.  The  particles  in  suspension  have  a  size 
distribution  and  their  relative  positions  are  described  by  a  pair-correlation 
function.  Twersky's  multiple  scattering  formalism  is  employed  to  obtain  new 
analytical  expressions  for  the  phase  velocity  and  coherent  attenuation  of  a 
wide  range  of  concentrations.  Numerical  results  presented  are  of  interest 
in  the  study  of  marine  sediments. 


INTRODUCTION 


Recent  advances  in  the  area  of  multiple  scattering  in  discrete  random 
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media  by  Twersky  have  shown  that  his  technique  can  uniformly  handle 
scattering  by  acoustic  and  electromagnetic  waves.  He  has  obtained  solutions 
in  the  long-wavelength  limit  for  the  bulk  propagation  constants  by  treating 
the  particles  as  "hard"  spheres  in  a  statistical-mechanical  sense.  In  this 
paper,  his  approach  is  applied  to  study  acoustic  wave  propagation  in  a  fluid 
containing  a  dense  suspension  of  spherical  elastic  particles  of  varying  sizes 
Closed  form  expressions  are  obtained  for  both  phase  velocity  and  coherent 
attenuation  in  the  Rayleigh  limit.  Numerical  results  are  presented  at  higher 
frequencies  for  marine  sediments  for  various  values  of  concentrations. 

SOLUTION  FOR  RADIALLY  SYMMETRIC  SCATTERERS 

Consider  a  random  distribution  of  elastic  spheres  in  water.  The  elastic 

properties  of  the  scatterers  are  given  by  Lame’s  constants  X^  and  p^  and 

density  p^.  The  fluid  properties  are  given  by  compressibility  X  and  density 

p.  The  longitudinal  and  shear  wave  velocities  in  the  elastic  scatterer  are 

1/2  1/2 

given  by  c  =  [ (X^  4-  2p^)/p^]  and  cg  =  (p^/p^)  while  the  compressional 

1/2 

wave  velocity  in  water  is  given  by  c  =  (X/p) 

The  effective  wavenumber  K  for  a  random  distribution  of  pair-correlated 

scatterers  can  be  written  in  terms  of  wavenumber  k(=  w/c)  in  water  as  follows 
Aim  00 

K2  =  k2  +  -j— ^  l  A  n2n  (1) 

k  L  n 
n=o 

where 

n  -  K/k 
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fn(a,0) 


fn(a,°)  q(a)da  ;  fn(a,0)  =  Y+Tc*1^ 

o  n 


Cn  is  a  complex  function  containing  density,  compressibility  and  rigidity  of 
the  elastic  scatterer  and  different  orders  of  spherical  Bessel  and  Neumann 
functions,  see  Appendix  A  and 


nm 
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In  Eqs.  (2)  and  (3),  nQ  =  N/V  is  the  number  density  (number  of  scatterers  per 
unit  volume),  fn(a,0)  is  the  forward  scattering  amplitude  of  a  spherical  scat¬ 
terer  of  radius  'a',  q(a)  is  the  size  distribution  function,  g(R)  is  the  pair- 

(2) 

correlation  functions,  )  and  hm  (  )  are  the  spherical  Bessel  and  Hankel 


functions,  respectively,  and  the  coefficients  d 


m 


(:■:) 


can  be  obtained  from 


the  expansion  of  the  Legendre  polynomials  as  given  by 


P  (x)  P  (x)  = 
n  m 


Pq(x)  . 


has  the  following  properties: 
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m+n  =  even 


=  1  ;  q  =  m+n 
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m+n  =  odd 


The  effective  wavenumber  K(=  K^+i^)  is  complex,  the  real  and  imaginary  parts 
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of  which  are  related  to  the  phase  velocity  and  coherent  attenuation,  res¬ 
pectively. 

For  monopole  and  dipole  cases 


In  the  Rayleigh  limit,  only  monopole  and  dipole  terms  and  a  two-term  expansion 
of  (1)  are  considered  in  obtaining  a  closed  form  expression  for  the  phase  veloc¬ 
ity  and  attenuation.  For  higher  frequencies,  the  general  form  of  d^  and  more 
terms  in  expansion  of  Eq.  (1)  should  be  employed.  This  is  best  done  numerically. 
To  calculate  the  effective  wavenumber  K,  the  equation  (3)  containing  N 

m 

which  involves  the  radial  distribution  function  g(R)  needs  to  be  obtained  first. 
For  a  more  precise  calculation,  the  correlation  integral,  Eq.  (3),  should  be 
computed  for  each  value  of  3.  However,  we  have  chosen  to  compute  the  corre¬ 
lation  integral  for  the  mean  sphere.  This  is  discussed  in  more  detail  in 

what  follows.  Following  the  work  by  TwersUy,  we  can  write 

r 

N  =  47m  F(K,k)ia  ;  ia  =  [g(R)-l]RadR  (7) 

mo 

o 

where  F(K,k)  is  a  function  of  the  effective  wavenumber  K  and  the  host  medium 
wavenumber  k.  Using  the  Virial  expansion'*'  for  g(R)  ,  we  can  show  that 


Ia  =  -(l+8w)  +  8w 


3)a+1  ,  3ba+1(l-2a+2)  ,  ba+1(2a+4-l) 
a+1  4(ct+2)  16(a+4) 


where  b  =  2a  and  w  =  n  4Tra/3  is  the  concentration  of  scatterers  and 

o 

a  -  \i/  4ttN/ 3w  is  the  mean  radius.  For  the  lowest  order  of  m  in  (7)  which  is 


useful  for  low  frequency  scattering,  N  is  obtained  as 

o 

=  47rn  1^  =  -  8w  +  34w^ 
o  .  o 


The  Virial  series  expansion  and  the  resulting  Eq.  (9)  seem  to  be  good  for  low 


concentrations  w<0.05.  For  higher  concentrations,  one  can  obtain  the  following 
expression  derived  by  Twersky  for  by  using  Ornstein-Zernicke  equation 


given  by 


N1  =  (1~  W-}  ■-  -  1  (1( 

°  (l+2w)2 

Eq.  (10)  is  preferred  to  Eq.  (9)  since  it  is  not  restricted  to  small  values  of 
w.  But  corresponding  expressions  cannot  be  obtained  for  all  values  of  a, 
hence  we  have  to  resort  to  the  virial  expansion. 


PHASE  VELOCITY  AND  COHERENT  ATTENUATION 


The  suspended  sedimentary  particles  in  water  are  in  general  2  to  4  microns 
in  radius.  For  an  acoustic  wave  with  a  frequency  range  of  1.2  to  12  MHz,  this 
corresponds  to  a  non-dimensional  wavenumber  ka  in  the  range  of  0.01  to  0.1. 

For  such  long  wavelength  acoustic  wave  propagation  studies,  it  is  enough  to 
keep  only  monopole  and  dipole  terms  in  Eq.  (1).  The  resulting  dispersion  equa¬ 
tion  is  given  by 
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Following  Twersky's  work'*',  Aq  and  A^  can  be  written  as  follows 
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For  ka"  <<  1,  we  obtain  the  following  limiting  values  for  aQ  and  a^ 
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In  Eqs.  (15)  and  (16),  g  is  the  density  ratio  between  the  scatterer  and  the 
surrounding  medium  and  c^,  cg  and  c  are  the  longitudinal  and  transverse  wave 
speeds  of  the  elastic  scatterer  and  the  compressional  wave  speed  of  the  sur¬ 
rounding  fluid.  As  a  limiting  case,  we  can  also  obtain  the  dispersion  rela¬ 
tions  for  a  random  distribution  of  either  rigid  (A^,Pj-*»)  or  fluid  (y^->0)  scat 
terers  in  water. 

After  some  manipulations,  one  can  obtain  the  real  (K^)  and  imaginary 


part  (K0)  of  the  effective  wavenumber  K  as  follows 
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In  Eq.  (17),  y  =  (l-8w+34w^)for  0<w<0.05  and  y= (1-w) ^/(l+2w) ^  for  w>0.0o. 
a  uniform  size  distribution,  G(aO=  1  while  for  a  Gaussian  size  distribution 

+  |(2  -  e*p(-  ^f)Kj)2 

+  —  exp[-  kf)  J(|)  +  kl  +  erf(-J->] 

/U  2  m  a  2  ^ 

where  m  is  the  standard  deviation. 

The  Gaussian  distribution  function  has  an  infinite  tail,  however  in 
the  numerical  calculations  the  cut  off  size  was  chosen  to  be  2a  where  a 
is  given  by 

a  =  ^/4ttN/3w 

With  this  cut  off,  for  ”  ■  0.4,  q(2a)  *  4.4x10  ^  which  is  small  (4%  of  peak 

a 

value) and  if  ka  is  small  2ka  is  still  quite  small  and  the  long  wavelength 


approximation  is  still  valid  for  m/a=0.18,  q(2a)=1.9xl0  ^ ,  which  is  negligible. 
The  question  may  be  raised  as  to  the  radius  of  the  excluded  spherical  volume 
surrounding  each  scatterer.  If  we  use  4a,  the  diameter  of  the  largest  size 
sphere  of  the  distribution,  although  we  would  be  correct  in  the  integration 
on  the  allowed  position  for  each  scatterer,  we  would  be  excluding  an  un¬ 
reasonable  volume,  thus  limiting  the  volume  fraction  of  scatterers.  It  seems 
more  reasonable  to  take  the  diameter  of  the  mean  sphere,  2a  as  the  radius 
of  the  excluded  spherical  volume.  On  the  average,  this  would  be  applicable 
to  most  scatterers  and  allow  us  to  consider  higher  volume  fractions.  If  the 
volume  fraction  is  low  enough,  whether  the  radius  of  the  excluded  sphere  is 
2a  or  4a  will  not  pose  a  problem.  For  higher  concentrations,  if  m/a  <0.5, 
the  procedure  we  have  followed  should  yield  good  results. 

Some  numerical  calculations  were  carried  out  for  sediment  particles  in 
water.  The  particles  are  assumed  to  be  spherical  in  shape  and  their  elastic 
properties  are  given  in  Table  I.  In  the  Rayleigh  limit,  the  normalized  phase 
velocity  K^/k  is  frequency  independent.  Its  dependence  on  concentration  is 
shown  in  Fig.  1  for  both  uniform  and  Gaussian  size  distributions.  The  co¬ 
herent  attenuation  (loss  tangent)  ^/K^,  however,  depends  on  frequency  as 
depicted  by  Fig.  2.  From  the  log-log  plot  of  Fig.  2,  one  can  see  that  the 
loss  tangent  for  both  uniform  and  gaussian  size  distributions  is  a  straight 
line  with  a  slope  equal  to  3  as  predicted  by  Eq .  (17).  Figure  3  presents 
the  loss  tangent  as  a  function  of  concentration  at  ka=0.1.  The  attenuation 
increases  with  increasing  concentration  until  it  reaches  a  maximum  at  c=15% 

and  then  decreases  with  increasing  concentration.  This  kind  of  behavior  has 

4 

been  observed  experimentally  by  Lloyd  and  Berry  . 

For  higher  frequencies  and  concentrations,  more  terms  must  be  kept  in 
Eq.  (1)  along  with  suitable  radial  distribution  function  g(R).  Our  model  of 
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the  random  system  is  that  the  spherical  scatterers  cannot  interpenetrate. 

In  the  statistical  mechanics  literature,  this  is  synonymous  with  the  ensemble 
of  "hard"  spheres.  Several  theories  and  calculations  are  available  for  deter¬ 
mining  the  joint  probability  distribution  function,  viz.,  the  Hypernetted-Chain 
Equation(HNC) ,  the  Percus-Yevick  Approximation  (P-YA) ,  the  Self-Consistent 
Approximation  (SCA),  Monte  Carlo  calculations,  etc.  We  found  the  SCA  is  better 
suited  for  higher  concentrations,  see  Ref.  5. 

In  Fig.  4,  the  normalized  phase  velocity  is  plotted  as  a  function  of 
frequency.  The  phase  velocity  displays  an  anomalous  dispersion  for  0.7<ka<1.0. 
This  anomaly  is  a  result  of  some  resonance  effect  associated  with  the  individual 
particles.  In  Fig.  5,  the  attenuation  is  plotted  as  a  function  of  frequency 
for  the  same  concentration,  c=0.105.  The  attenuation  is  a  maximum  in  the  fre¬ 
quency  range  where  the  phase  velocity  displays  an  anomalous  behavior.  However, 
since  the  calculations  were  made  only  up  to  ka=1.0,  we  cannot  make  a  definitive 
statement  in  this  regard.  In  Figs.  6  and  7,  the  normalized  phase  velocity  and 
attenuation  are  plotted  as  a  function  of  concentration  at  ka=1.0. 


CONCLUSIONS 


In  this  paper,  Twersky*s  multiple  scattering  formalism  has  been  ex¬ 
tended  to  elastic  scatterers  in  water.  It  has  been  shown  that  in  addition 
to  providing  closed  form  expressions  for  the  dispersion  relations  for  elastic 
scatterers  in  the  low  frequency  limit,  his  formalism  is  also  well  suited  for 
numerical  computations  at  higher  frequencies  and  concentrations  employing 
suitable  pair-correlation  functions. 
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TABLE  I. 


Elastic  Properties  of  Sediments  Used  in  Calculations 
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Range  of  grain  sizes  (a) 

Grain  density  (p^) 

Shear  wave  velocity  (cg) 
Compressional  wave  velocity  (c  ) 


2-4  pm 
2.65  g/cm^ 
210  m/sec 
3690  m/sec 


Sound  speed  in  water 


1500  m/sec 
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APPENDIX  A:  THE  EXPRESSION  OF  C  USED  IN  EQ.  (1) 
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APPENDIX  A  (Cont.) 
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nn  =the  n-th  order  spherical  Newman  function,  j  =the  n-th  order  spherical 

Bessel  function. 
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FIGURE  CAPTIONS 


Figure  1.  Phase  velocity  vs  concentration  w  in  the  Rayleigh  limit  for 
Gaussian  (m/a  =  0.18,  0.4)  and  uniform  size  (m/a  =  0.0)  dis¬ 
tributions 

Figure  2.  Coherent  wave  attenuation  vs  ka  for  w=0.1  and  m/a  =  0.0,  0.18,  0.4 

Figure  3.  Coherent  wave  attenuation  vs  concentration  w  for  ka  =  0.1  and 
m/a  =  0.0,  0.18,  0.4. 


Figure  4. 
Figure  5. 
Figure  6. 
Figure  7. 


Normalized  phase  velocity  vs  ka  for  w=0.105. 
Coherent  wave  attenuation  vs  ka  for  w=0.105. 
Normalized  phase  velocity  vs  w  for  ka  =  1.0. 
Coherent  wave  attenuation  vs  w  for  ka  =  1.0. 
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for  very  high  frequencies *(ka>50)  r~~-The  theoretical  results  obtained  by  this 
formalism  were  compared  with  experimental  findings  and  the  agreement  is 
excellent. 
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